IBRARY 


■AVAL  POS'iGRAOUAH  SCHOOL 

lONTERtY.  CALIFORNIA  93940 


yWRDC-TR-8 9-3089 


(FREQUENCY  DOMAIN  DESIGN  OF  ROBUST  CONTROLLERS 

'For  space  structures 


AM  ZH 


W.H.  Bennett 
C.  LaVigna 

i_e  \  Systems  Engineering,  Inc.  >Tt IT  -  |  S' -  V/Jfi 

7833  Walker  Dr.,  Suite  620 
Greenbelt,  MD  20770 


l! 


August 


1989 


Final  Report  for  Period  September  1988  -  February  1989 


Approved  for  Public  Release;  Distribution  is  Unlimited 


FLIGHT  DYNAMICS  LABORATORY 
WRIGHT  RESEARCH  AND  DEVELOPMENT  CENTER  - 
'fAIR  FORCE  SYSTEMS  COMMAND 
WRIGHT-PATTERSON  AIR  FORCE  BASE,  OHIO  45433-6553 


NOTICE 


When  Government  drawings,  specifications,  or  other  data  are  used  for 
any  purpose  other  than  in  connection  with  a  definitely  Government-related 
procurement,  the  United  States  Government  incurs  no  responsibility  or  any 
obligation  whatsoever.  The  fact  that  the  government  may  have  formulated  or 
in  any  way  supplied  the  said  drawings,  specifications,  or  other  data,  is  not 
to  be  regarded  by  implication,  or  otherwise  in  any  manner  construed,  as 
licensing  the  holder,  or  any  other  person  or  corporation;  or  as  conveying 
any  rights  or  permission  to  manufacture,  use,  or  sell  any  patented  invention 
that  may  in  any  way  be  related  thereto. 


This  report  is  releasable  to  the  National  Technical  Information  Service 
(NTIS).  At  NTIS,  it  will  be  available  to  the  general  public,  including 
foreign  nations. 


This  technical  report  has  been  reviewed  and  is  approved  for  publica 

tion . 


Projeci  Engineer 

Design  A'  Analysis  Methods  Group 
FOR  THE  COMMANDER 


NELSON  D.  WOLF,  Technical  Manager 
Design  k  Analysis  Methods  Group 
Analysis  k  Optimization  Branch 


_rOHN  T.  ACH,  Chief 
Analysis  Optimization  Branch 
Structures  Division 


If  your  address  has  changed,  if  you  wish  to  be  removed  from  our  mailing 
list,  or  if  the  addressee  is  no  longer  employed  by  your  organization  please 
notify  WRDC/FIBRA.  WPAFB,  OH  45433-  6553  to  help  us  maintain  a  current 
mailing  list. 


Copies  of  this  report  should  not  be  returned  unless  return  is  required  by 
security  considerations,  contractual  obligations,  or  notice  on  a  specific 
document . 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  THIS  PAGE 


REPORT  DOCUMENTATION  PAGE 


la  REPORT  SECURITY  CLASSIFICATION 

UNCLASSIFIED 


lb.  RESTRICTIVE  MARKINGS 

NONE 


2a.  SECURITY  CLASSIFICATION  AUTHORITY 


2b.  DEC  LA  SSI  F  I C  A  T I  ON/DOWN  GRADING  SCHEDULE 


3.  DISTRIBUTION/AVAILABILITY  OF  REPORT 

Approved  for  public  release;  distribution 
is  unlimited. 


4  PERFORMING  ORGANIZATION  REPORT  NUMBE  R  (S) 

SEI-89-03-15-WB 


5.  MONITORING  ORGANIZATION  REPORT  NUMBER(S) 

WRDC-TR-8 9-3089 


6a.  NAME  OF  PERFORMING  ORGANIZATION 

Systems  Engineering,  Inc. 


6b.  OFFICE  SYMBOL 
(If  applicable ) 


7a.  NAME  OF  MONITORING  ORGANIZATION 

Air  Force  Wright  Research  &  Development 
Center  (WRDC/F1BRA) 


6c.  ADDRESS  (City,  State  and  ZIP  Code) 

7833  Walker  Dr.  Suite  620 
Greenbe It  MD  20770 


7b.  ADDRESS  (City.  State  and  ZIP  Code) 


Wright-Patterson  Air  Force  Base,  Ohio 

45433-6553 


8a  NAME  OF  FUNDING/SPONSORING 
ORGANIZATION 

Strategic  Defense  Initiative 


Bb.  OFFICE  SYMBOL 
(If  appticable) 


9.  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 

F33615-88-C-3215 


8c.  ADDRESS  (City.  State  and  ZIP  Code) 

T/IS/SBIR 


10.  SOURCE  OF  FUNDING  NOS. 


11  TITLE  (Include  Security  Classification ) 

Frequency  Domain  Design  of  Robust  Controllers 


PROGRAM 
ELEMENT  NO 

63220C 

for  Space  St 


PROJECT 

TASK 

WORK  UNIT 

NO 

NO 

NO. 

2200 

51 

89 

ructures 

12.  PERSONAL  AUTHOR (S) 

W.H.  Bennett  and  C.  LaVigna 


13a.  TYPE  OF  REPORT 

13b.  TIME  COVERED 

14  DATE  OF 

REPORT  (Yr  .  Mo.,  Day) 

15.  PAGE  COUNT 

FINAL 

from  Sep  88  to  Feb  89 

15 

89 

53 

16.  SUPPLEMENTARY  NOTATION 


17  COSATI  CODES 

FIELD 

GROUP 

SUB.  GR. 

IB.  SUBJECT  TERMS  (Continue  on  reverse  if  necessary  and  identify  by  Mock  number) 

SEE  ATTACHED 


19.  ABSTRACT  (Continue  on  reverse  if  necessary  and  identify  by  block  number ) 


SEE  ATTACHED 


20.  DISTRIBUTION/AVAILABILITY  OF  ABSTRACT 
UNCLASSIFIED/UNLIMITED  3  SAME  AS  RPT.  □  OTIC  USERS  □ 


21  ABSTRACT  SECURITY  CLASSIFICATION 

UNCLASSIFIED 


22a.  NAME  OF  RESPONSIBLE  INDIVIDUAL 

Dr.  V.  Venkayya 


22b.  TELEPHONE  NUMBER 
(Include  Area  Code) 

(513)255-7191 


22c.  OFFICE  SYMBOL 

WRDC/FIBRA 


DD  FORM  1473,  83  APR 


EDITION  OF  1  JAN  73  IS  OBSOLETE. 


-UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  THIS  PAGE 


Small  Business  Innovation  Research  (SBIR)  Program 
Phase  1  -  Summary  Report 

Contract  No.  F33615-88-C-3215 


Topic  No.:  SDI87-12 

Small  Business  Firm:  SYSTEMS  ENGINEERING,  INC.1 
Principal  Investigator:  Dr.  William  H.  Bennett 

Proposal  Title:  Frequency  Domain  Design  of  Robust  Controllers  for  Space 
Structures 


Project  Summary: 

The  purpose  of  the  Phase  1  effort  was  to  investigate  and  demonstrate  the 
feasibility  of  a  new  class  of  computational  algorithms  for  the  design  of  high 
performance  control  laws  for  flexible  space  structures  based  on  frequency 
response  modeling  and  to  consider  advanced  techniques  for  the  implementa¬ 
tion  of  real  time  control  for  precision  (wide  bandwidth)  applications.  Typical 
applications  requiring  advanced  realtime  control  of  flexible  space  structure 
include  vibration  suppression  and  isolation  of  payload  subsystems.  Perfor¬ 
mance  of  vibration  suppression  and  isolation  systems  are  critical  factors  ef¬ 
fecting  achievable  levels  of  performance  for  space  based  optical  systems. 

In  the  Phase  1  effort,  a  prototype  software  code  was  developed  for  testing 
computational  algorithms  for  spectral  factorization,  causal  projection,  and 
coprime  factorization — critical  steps  in  frequency  domain  design  of  precision 
control  laws.  Several  representative  models  of  vibration  suppression  and  iso¬ 
lation  of  flexible  structures  were  developed  and  control  laws  were  designed 
and  tested.  The  innovative  approach  for  frequency  domain  computations 
employed  is  based  on  sampling  and  interpolation  of  the  system  frequency 
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numerical  stability  of  the  computational  algorithms  for  both  irrational  (dis¬ 
tributed  parameter)  transfer  function  models  and  various  rational  (lumped 
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1  Phase  1  Project  Objectives  and  Progress  Summary 

Available  methods  for  design  of  complex,  multiloop  control  systems  for  flexible  structures  are 
currently  quite  limited.  Current  methods  are  largely  based  on  state  space  models  obtained  by 
modal  truncation  and  controllers  are  designed  by  LQG  optimization  [l].  Achieving  precision 
control  will  depend  on  the  ability  to  optimally  assess  design  tradeoffs  for  models  which 
capture  the  essential  distributed  parameter  nature  of  the  dynamic  phenomenon  involved  [2]. 

The  Modern  Control  Theory  Paradigm:  Computation  and  Implementation  by 
State  Space  Methods  Modern  methods  in  control  system  design  are  based  on  a  class  of 
mathematical  models  for  the  linear  (small  amplitude)  response  of  systems  (called  state  space 
models)  resulting  from  modeling  the  system  response  as  a  system  on  first  order  Ordinary 
Differential  Equations  (ODE’s)  in  time.  Optimization  of  the  system  time  responses  is  used  to 
resolve  certain  basic  engineering  tradeoffs.  State  space  models  are  particularly  significant  for 
time  response  optimization  since  they  provide  an  internal  realization  of  the  system  dynamics 
in  terms  of  a  system  state  which  codifies  the  past  history  of  the  system  response  up  to  the 
current  time.  The  concept  of  a  state  space  realization  has  proven  advantageous  (at  least 
for  low-order  systems)  for  several  reasons  including  the  fact  the  computational  algorithms 
for  time  response  optimization  can  be  described  in  terms  of  constant  real  matrices.  Such 
computational  algorithms  can  be  readily  supported  by  realiable  numerical  software  running 
on  a  standard  digital  computer. 

In  the  theory  the  system  response  is  modeled  by  a  state  space  model; 

x(t)  =  Ax(t)  +  Bu(t)y 
y(t)  =  Cx(t), 

where  x(t)  £  3Rn  is  the  system  state,  y(£)  £  is  the  vectyor  of  system  outputs,  and 
u(£)  £  is  the  vector  of  system  inputs.  Control  system  design  methods  are  based  on 
the  separation  of  state  variable  feedback  and  asymptotic  state  estimator  (observer)  designs. 
The  resulting  control  compensator  suggested  by  this  process  of  separate  design  steps  can  be 
realized  as  shown  in  Figure  1.1.  Here  the  triangular  blocks  designate  an  integration  operation 
for  continuous  time  system  realizations  as  described  above.  The  picture  (drawn  in  this  way) 
is  suggestive  of  the  principal  technology  for  simulation  of  dynamical  systems — popular  in 
the  1950-  1970’s — associated  with  analog  or  hybrid  computers. 

As  the  performance  of  digital  computers  increased  while  their  size  and  power  require¬ 
ments  decreased  a  range  of  applications  for  realtime  control  using  digital  computers  became 
possible.  The  modern  control  paradigm  was  then  extended  by  considering  the  state  space 
realization  for  the  equivalent  sampled-data  system  of  the  form 

x{(k  +  1)T)  =  Ax{kT)  +  Bu(kT), 
y{kT)  =  Cx(fcT), 

where  T  is  the  (uniform)  sampling  interval.  In  this  case  the  implemenation/realization 
paradigm  carries  over  directly  by  considering  the  triangular  blocks  in  Fig.  1.1  to  be  a  unit 
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Figure  1.1:  The  realization/implementation  paradigm  of  modern  control  theory. 

delay  of  time  T.  However,  it  is  apparent  that  such  a  component  based  approach  to  imple¬ 
mentation  offers  a  constrained  view  of  options  for  realtime  control  implemenation  in  a  digital 
computer. 

It  is  the  basic  proposition  of  this  project  (and  the  phase  2  proposal  will  focus  on  this 
issue)  that  the  above  component  based  approach  to  control  law  implementation  by  state 
space  realization  is  essentially  obsolete  given  the  state-of-the-art  in  high  speed  VLSI  designed 
for  signal  processing  (i.e.  read-time  applications.)  We  contend  that  the  various  technologies 
for  high  speed  and  very  wide  time-bandwidth  product  computations  (including  advanced 
VLSI  and  special  acousto-optical  systems)  are  now  mature  enough  to  investigate  alternate 
paradigms  for  readtime  control  law  implementation. 

State  Space  Alternatives  for  Distributed  Parameter  Systems.  For  application  to 
control  design  computations  for  flexible  space  structures  we  note  the  use  of  state  space  models 
may  be  particularly  awkward  for  the  typical  computations  required  to  obtain  time  responses. 
This  is  true  in  general  for  distributed  parameter  systems  since  the  appropriate  state  space 
models  will  involve  a  coupled  set  of  integral-partial  differential  equations  (IPDE’s).  However, 
for  the  class  of  models  typically  encountered  in  dynamics  of  space  structures  we  observe  that 
transfer  function  (or  frequency  response)  models  may  be  particularly  significant. 

We  suggest  that  the  use  of  transfer  function  models  for  control  design  can  offer  significant 
alternatives  for  control  system  design  for  the  following  reasons. 

1.  The  interplay  between  ‘state-space’  and  ‘input-output’  methods — evident  in  the  theory 
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of  finite-dimensional  systems  [3,  4] — is  even  more  important  for  distributed  parameter 
systems  with  potential  for  profound  results  in  basic  engineering  methods. 

2.  Precise  engineering  specifications  for  distributed  control  system  designs  are  extremely 
difficult  to  formulate  in  terms  of  state-space  models  (e.g.,  in  terms  of  PDE’s  and  their 
coefficients).  These  specifications  can  be  rather  easily  characterized  in  terms  of  the 
system  frequency  response. 

3.  By  working  with  the  transfer  function  model  directly,  alternate  control  law  implemen¬ 
tations  can  be  discovered  which  may  be  simple  to  implement  and  very  effective. 

4.  It  is  convenient  to  characterize  model  uncertainty  in  a  frequency  domain  setting  and 
this  approach  provides  insight  useful  for  control  system  design  [5]. 

Remarks  on  Implementation  of  Distributed  Control  For  a  variety  of  reasons  dis¬ 
tributed  state  feedback  control  laws  may  be  difficult  to  implement  and  state  observer-based 
control  will  require  the  on-line  realization  of  a  distributed  parameter  system.  In  the  mod¬ 
ern,  state-space  approach  to  design  the  control  law  realization  is  linked  with  the  state  space 
model  for  the  system.  We  assert  that  this  paradigm  may  seriously  restrict  design  options 
for  control  law  realization  for  a  large  class  of  linear,  time-invariant  systems  including  dis¬ 
tributed  parameter  effects.  We  prefer  to  think  more  generally  of  a  linear,  dynamic  control 
law  as  modeled  by  a  convolution  equation, 

u(t)  =  C(t)*y(t), 

where  y(t)  are  available  sensor  measurements  and  u(t )  are  control  actuation  signals.  For 
linear,  time-invariant  system  models  the  control  law  can  be  specified  by  its  impulse  response 
0(t)  or  in  the  frequency  domain  as  a  transfer  function  (7(.s)  which  may  be  irrational.  The 
principal  innovation  of  this  project  is  that  control  system  design  for  distributed  parame¬ 
ter  systems  may  benefit  from  decoupling  the  computation  and  realization  design  steps.  A 
principal  project  objective  of  this  study  is  to  assess  the  capability  for  alternate  methods  for 
control  law  realization  which  may  offer  reasonable  approximation  of  an  optimal  irrational 
transfer  function  for  a  control  law. 


2  Technical  Background  for  Project  Innovation 


2.1  Transfer  Function  Models  for  Flexible  Structures 


Generic  models  for  the  elastic  response  of  flexible  structures  are  often  described  as  spatial 
continuum  via  a  Partial  Differential  Equation  (PDE)  of  the  form, 


m(z ) 


d2w(t,z)  t  n  diu(t,  z) 

dt*  0  at 


+  A0u>(t,  z)  =  F(t,  z) 


(2.1) 


where  to(t,z)  is  an  N-vector  of  displacements  of  a  structure  17  with  respect  to  some  equi¬ 
librium  for  ft  is  a  bounded,  open  set  in  9?^  [1].  The  (vector)  z  €  17  is  a  coordinate  in  17. 
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We  assume  the  boundary  <9f 2  is  smooth.  The  mass  density  m(z)  is  positive  definite  and 
bounded  on  dQ.  The  damping  term  Dodw/dt  models  both  (asymmetric)  gyroscopic  and 
(symmetric)  structural  damping  effects.  The  internal  restoring  force  A0iu  is  generated  by 
a  time-invariant,  differential  operator  A0  for  the  structure.  For  most  common  structural 
models,  A0  is  an  unbounded,  differential  operator  with  domain  D(A0)  consisting  of  certain 
smooth  functions  satisfying  appropriate  boundary  conditions  on  <9fi.  Thus,  for  these  prob¬ 
lems,  D(A0)  is  typically  dense  in  the  Hilbert  space  Ho  =  £2^)  endowed  with  its  natural 
inner  product  (x,t/)0  =  /q  xT(z)y(z)  dz.  Often  (but  not  always),  the  spectrum  of  A0 ,  <j(Aq), 
consists  of  discrete  eigenvalues  with  associated  eigenfunctions  which  constitute  a  basis  for 

c2(Q). 

The  applied  force  distribution  F(£,  z)  can  be  thought  of  as  consisting  of  three  components 

F(ty  z)  =  Fd(t,  z)  +  Fc(t ,  z )  +  Fa(t ,  z)  (2.2) 

where  F ^  is  N- vector  of  exogenous  disturbances  (possibly  forces  and  torques),  Fc  is  a  continu¬ 
ous,  distributed  controlled  force  field  (an  available  option  in  only  some  special  applications), 
and  Fa  represents  controlled  forces  due  to  localized  actuation; 

k 

Fa(t ,  z)  -  ^2  bj(z)uj(t)  =  B0u(t).  (2.3) 

j=t 

The  actuator  influence  functions  bj(z)  are  highly  localized  in  and  can  be  approximated  by 
Dirac  delta  functions.  We  assume  that  a  finite  number  p  of  measurements  can  be  made  as: 

j,(f)  =  Cow  +  C'0jj-  (2.4) 

where  y(£)  is  a  p- vector.  The  operators  B0  :  7 i0,  Co  •  FLq  — > ►  and  Cq  :  Ho  — ►  5?p 

are  bounded. 

A  natural  assumption  for  structural  problems  [1]  is  that  A0  is  self-adjoint  with  compact 

resolvent  and  discrete  (real)  spectrum  which  is  bounded  from  below.  The  system  state  can 

1 

then  be  considered  as  an  element  of  a  Hilbert  space  H  =  D(Aq  )  x  Ho  with  the  energy  norm 

IMIl  =  («b  -^0^)0  +  {miv,w)0  (2.5) 

where  the  first  term  represents  potential  energy  and  the  second  term  is  kinetic  energy.  Thus 
the  (abstract)  state  space  model  can  be  written 

x(tyz)  =  Ax(t,  z)  +  Bu{t)  (2.6) 

y{t)  =  Cx(t,z) 


where 


I 

-Do 


c  =  [c„,  c;i . 


(2.7) 


For  the  elastic  dynamics  of  space  structures,  there  is  always  some  (possibly  small)  damping 
Dq  appearing  in  (2.1)  which  causes  A  to  be  dissipative.  Thus  the  criteria  of  the  Hille-Yoshida- 
Phillips  theorem  [6]  are  satisfied  and  A  generates  a  (70-seniigroup  with  an  operator  which 
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we  write  suggestively  as  eAt.  Moreover,  such  models  are  ‘hyperbolic’  [1]  in  the  sense  that 
the  semigroup  is  a  contraction,  i.e.,  || ||  <  1  and  all  but  the  zero  frequency  poles  are  only 
slightly  damped;  ||ej4t||  <  e~St  for  some  small  5  >  0.  We  remark  that  some  popular  models 
for  structural  elements  such  as  beams  with  material  damping  may  not  fit  in  this  framework 
[7,  8].  However,  this  framework  includes  models  appropriate  for  considerations  of  wave-like 
dynamics  which  propagate  causally  in  the  spatial  domain.  For  such  models,  the  question  of 
how  to  compute  controls  u(£)  and  system  response  x(£)  focuses  on  the  so-called  weak  solution 
of  (2.6); 

x(£,  z)  =  eAtx(0,  z)  +  f  Bu(cr)d<T.  (2.8) 

Jo 

PDE’s  encountered  in  generic  models  of  the  form  (2.1)  represent  the  time  evolution  of 
certain  physical  systems  and  are  typically  of  either  hyperbolic  or  parabolic  type.  To  clarify 
terminology  consider  the  system  of  first-order,  partial  differential  equations  defined  for  t  >  0 
and  0  <  z  <  L 

dw  -  dw 

ETt=FTz  +  Hw '  w  s  fl"  <29> 

If  E  is  nonsingular,  then  (2.9)  can  be  written 


dw  dw 
~dt  =  F~d~z 


+  Hw 


(2.10) 


where  F  =  E~1F1  H  —  E~lH.  If  F  has  only  real  eigenvalues  and  a  complete  set  of 
eigenvectors,  then  the  system  is  said  to  be  hyperbolic  (see,  for  example,  Zauderer  (9]).  If 
there  are  multiple  real  eigenvalues  and  less  than  a  complete  set  of  eigenvectors,  then  the 
system  is  of  (partial)  parabolic  type.  If  all  of  the  eigenvalues  are  complex,  the  system  is  of 
elliptic  type.  Systems  with  complex  eigenvalues  are  not  causal  and  will  not  be  considered 
further. 

If  E  is  singular,  (2.9)  can  give  rise  to  mixed  systems  of  all  types.  Our  interest  in  this 
case  will  be  limited  to  purely  parabolic  systems  of  the  type 


dw  rid2w 
——  =  Gr  * 


_dw 

+  F—  +  Hw 
oz 


(2,11) 


dt  ~  dz1 

which  commonly  arise  in  engineering  problems. 

In  addition  to  equations  (2.9)  or  (2.11),  there  are  associated  initial  and  boundary  condi¬ 
tions.  For  equation  (2.9),  these  conditions  take  the  general  form 


initial  conditions  u>(z,0)  =  f(z) 

boundary  conditions  £iu>(0,£)  +  Tiw(L,t)  =  g(t) 


(2.12) 


where  dim(y)=dim(tw);  and  for  equation  (2.11),  they  take  the  general  form 
initial  conditions  ui(z,0)  =  /(z) 

boundary  conditions  Eitt>(0 ,t)  +  E2|^(0,t)  +  T\w(L,t)  +  T2^z(L,t)  =  g(t) 
where  dim(^)=2dim(w). 

It  is  well  known  that  the  coefficient  matrices  in  (2.12),  (2.13)  must  satisfy  certain  con¬ 
straints  if  the  problem  formulation  is  to  be  well-posed.  In  the  hyperbolic  case  (equations  (2.9) 
and  (2.12)),  these  constraints  essentially  require  that  the  boundary  conditions  be  compatible 
with  the  wave  directions. 
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2.1.1  Frequency  Response  Models  for  Elastic  Structures 

The  class  of  models  discussed  above  captures  the  time  response  of  causal,  linear  systems 
arising  in  the  modeling  of  flexible  structures.  For  the  computation  of  optimal  control  laws 
for  such  systems  we  are  concerned  with  models  for  the  effective  system  frequency  response. 
A  complete  modeling  approach  is  described  in  [2]  which  we  summarize  in  the  next  few 
paragraphs. 

The  complete  system  transient  response  can  be  characterized  in  the  frequency  domain 
in  terms  of  the  superposition  of  the  forced  response  (modeled  via  a  transfer  function)  and 
the  free  response  (modeled  via  a  Green’s  function).  To  illustrate  the  approach  consider  a 
hyperbolic  model  in  one  space  dimension  0  <  z  <  L, 

tm=F?^f+Hx{t,z)+Ev{l,z)  (2.i4) 

subject  to  boundary  conditions 

£iz(*,0)  +  =  Df(t ),  (2.15) 

and  initial  conditions 

x{0,z)  =  x°(z)eHn{0,L).  (2.16) 

Here,  x  is  an  n-vector  valued  state  x  6  7in(0,L),  v  £  7il(0,  L)  is  an  ^-vector  valued  dis- 
tributed  disturbance,  /  is  m-vector  valued  boundary  interactions,  F,  H  are  real  n  x  n 
matrices  with  F  nonsingular  and  diagonalizable  [9],  and  Ej,  T\  are  n  x  n  matrices.  After 
taking  Laplace  transforms  in  the  temporal  variable  t,  we  obtain 

rL 

X($,z)=  /  Gr{s,  Z,  w)M(sjw)dw  +  Hbc{s,  z)F{s)  (2.17) 

Jo 

where 

M{s ,  z)  —  z°(u;)  —  C  u;), 

and  Xy  V ,  F  are  the  Laplace  transforms  of  x,  v,  /  respectively.  The  function  G>(s,z,tu) 
is  the  Greens  function  [10,  11]  for  (2.14),  (2.15)  and  Hbc{s)z)  is  a  transfer  function  from 
boundary  interactions  to  state.  Since  in  most  cases  of  practical  interest  the  control  of  flexible 
structures  will  be  eflected  by  actuators  whose  influence  functions  are  highly  localized,  we 
have  formulated  our  model  with  boundary  control  only. 

The  system  transfer  function  with  respect  to  boundary  control  has  the  form 

Hbc{s,z)  -  N(s,z)D  (2.18) 

where 


W(j,z)  =  #(j,z)|S,  +rlf(s,i)]-1,  (2.19) 

The  Green’s  function  for  (2.14),  (2.15)  is  the  solution  to 

— ^7T  — ~  =  F~l  [•»/  -  H ]  Gr(s,  z,  w)  +  I„S(z  -  w)  (2.20) 
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subject  to  the  boundary  conditions 


SiGr(s,  0,  w)  4-  TiGr(s,  L,w)  =  0 


(2.21) 


where  S(-)  is  the  Dirac  delta  function  [10],  [11].  Further  discussion  of  the  computation  of  the 
required  Green’s  functions  is  discussed  in  [2]. 

For  the  case  of  a  parabolic  model  we  begin  with  the  system  equations  in  the  form, 


dx(t,  z ) 

ft 


z) 


+  F 


dxjt, 

dz 


z) 


+  Hx(t ,  z)  +  Ev(t,  z) 


subject  to  2 n  boundary  conditions 


(2.22) 


Si*(t,0)  +  +  r 1x(tt  L)  +  =  Df(t),  (2.23) 

and  n  initial  conditions 

x{0,z)  =  X°(z)e  Hn{0,L).  (2.24) 

Let  £  =  [S1,S2]  and  V  =  [T i ,  T2] — each  2n  x  2n  real  matrics.  The  transfer  function  from 
boundary  control  is  then  of  the  form 


Hbc(3>z)  =  M(s,z)D, 

where 

*  /  \  _  On  In 

A\3>  =  -  sln)  -G-'F  j  ’ 

$(s,z)  =  eA^% 

M(s,z)  =  [/^Ol^.zMS  +  r^s,!)]-1. 


(2.25) 


(2.26) 

(2.27) 


Hybrid  Models  for  Flexible  Space  Structures.  In  most  applications,  models  for  the 
dynamics  of  flexible  structures  involve  interaction  between  various  elastic  and  rigid  compo¬ 
nents.  In  the  particular  case  of  flexible  structures  associated  with  large  space  structures,  the 
potential  topological  configurations  can  be  quite  complex.  Various  elements  such  as  beams, 
truss  structures,  cables,  membranes,  etc.,  may  have  dominant  distributed  parameter  effects. 
Typically  a  central  body  or  bodies  represent  large  concentrations  of  mass  with  respect  to 
the  overall  low  mass  density  of  the  flexible  structure.  These  are  most  effectively  represented 
by  lumped  parameter  models  of  their  rigid  body  dynamics.  Hybrid  models  can  provide  an 
effective  tool  for  analysis  of  dynamics  of  vibrations  and  their  effect  on  small  angle  motions 
for  complex  space  platforms. 

To  assemble  the  hybrid  model  we  consider  the  Distributed  Parameter  System  (DPS)  to 
be  modeled  as  either  the  hyperbolic  or  parabolic  (or  mixed)  cases  which,  as  we  have  seen, 
can  be  expressed  in  the  frequency  domain  in  the  form 

tL 

X d(s,z)=  /  Gr{siz1vj)M(31vj)dw  +  HBc{s,  q)Fd(s) 

Jo 


(2.28) 
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where 

M(s,u;)  =  x2(w)  “  EV(s,w).  (2.29) 

Clearly,  ( 2.28 )-( 2.29)  can  represent  a  disjoint  collection  of  distributed  elements  such  as  beams, 
cables,  etc. 

The  Lumped  Parameter  System  (LPS)  component  model  are  combined  into  a  single  LPS 
state  model  of  the  form 


*i{t)  =  Aixt(t)  +  x°t  =  Zi{0)  (2.30) 

with  x i  £  RN*  =  Xt  a  finite  dimensional  real  space.  By  taking  Laplace  transforms  in  (2.30), 
we  write  (analogous  to  (2.28)) 

Xtis)  =  Rt(s)z°t  +  Ht(s)F((s),  (2.31) 

where  R({s)  =  At]-1  is  the  resolvent  for  the  (matrix)  operator  At  and  Ht(s)  =  R((s)B. 

Finally,  the  interconnection  of  component  systems  is  resolved  through  a  topological  con¬ 
straint  relation  consisting  of  m  =  m.j  +  ra*  linear  equations; 

/(0  +  ri*rf(<,o)  +  r,*rf(4lL)  +  r,*<(0  =  Ku{t)  (2.32) 


where  u(t)  is  a  A:-vector  of  control  inputs  to  the  hybrid  system,  T\,T2  are  m  x  Nd ,  Ty  is 
m  x  N(,  and  K  is  m  x  k  real  matrices.  The  hybrid  model  state  space  consists  of  elements  of 
the  form 


x(t,z) 


\ 

Xd{t,z)  ) 


(2.33) 


which  are  N  =  Nd  +  TVr valued  functions  of  z  £  [0,  L]  ,  t  >  0.  The  resulting  model  has  the 
form 


rL  _  A  _  . 

A"(,s,z)  =  /  Gr(s}  z,w)M(s,w)dw  +  R(s,  z)x°t  +  Hbc{*,  *)&{*),  (2.34) 

Jo 


where  M(s,w)  is  given  in  (2.29).  The  complete  system  response  is  given  in  terms  of  the  free 
response  to  initial  conditions  in  the  LPS  state  via  an  N  x  matrix  transfer  function  and 
the  free  response  to  initial  conditions  in  the  DPS  state  given  by  a  N  X  Nd  Green’s  function 

(X; 


6/  Nc(s)Q \(s) 

R(,'z)  -  J  ™a)' 

(2.35) 

G,(,,z,w)  =  -#<(.)<?,(»)  1 />(,,„) 

[  Gr(3,  z,  IV)  -  Hbc{s,  z)Q2[s)  J 

(2.36) 

<?(»)  =  [/«  +  <?(<-)!-'  =  f  1 , 

.  V2(>»)  . 

(2.37) 

Q{*)  =  [TsH^^HBcis^  +  TtHBc^L)}, 

(2.38) 

P{s,iu)  =  TxGr(s,0,w) +  T2Gt(s,L,w). 

(2.39) 

where 
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The  forred  response  is  given  by  the  transfer  function  from  boundary  control, 


Ht(s)  0 
0  Hgc{s,z) 


Q(*)K. 


(2.40) 


In  the  phase  1  effort  we  considered  the  transfer  function  models  of  several  simple  sys¬ 
tems  respresentative  of  vibration  control  problems.  We  would  like  to  emphasize  that  the 
control  methods  we  have  considered  and  tested  in  Phase  1  are  by  no  means  limited  to  DPS 
applications.  In  fact  we  illustrate  the  computations  required  using  a  standard  finite  element 
model  in  one  of  the  benchmark  problems  considered.  The  significant  aspect  of  the  modeling 
formulation  above  is  in  the  fact  that  we  require  only  frequency  response  models  and  there¬ 
fore  the  computational  approach  is  general  enough  to  include  a  class  of  DPS  large  enough 
to  encompass  flexible  space  structure  applications. 


2.2  Control  Synthesis  via  Wiener-Hopf  Methods 

A  complete  frequency  domain  design  method  for  multi-input /multi-output  (MIMO)  systems 
based  on  a  optimal  Wiener-Hopf  problem  was  first  given  by  Youla  et  al  [3]  in  1976.  The 
method  described  by  Youla  and  his  colleagues  are  very  significant  for  control  system  design 
for  several  reasons.  The  method  directly  addresses  several  key  design  objectives  and  tradeoffs 
for  closed  loop  control  and  uses  the  setting  of  Wiener-Hopf  optimization  to  resolve  tradeoffs. 
The  formulation  of  the  problem  in  the  frequency  domain  permits  both  stabilization  and 
system  performance  to  be  addressed  in  the  same  objective.  This  is  accomplished  by  the 
characterization  of  an  algebraic  condition  depending  on  the  transfer  function  of  the  plant 
which  must  be  satisfied  by  any  closed  loop  controller  which  stabilizes  this  plant.  Then 
performance  objectives  can  be  defined  in  a  frequency  dependent  setting  without  regard  to 
stability.  The  solution  of  the  resulting  optimal  control  problem  is  obtained  by  solving  a  set  of 
Wiener-Hopf  problems  subject  to  the  constraint  that  the  controller  must  stabilize  the  plant. 

Youla’s  approach  focuses  on  a  specific  control  architecture  displayed  in  Figure  2.1  of  prac¬ 
tical  significance  in  a  wide  variety  of  control  problems  where  P(s)  is  an  n  x  m  plant  transfer 
function  and  F(s)  is  n  X  n  and  models  sensor  measurement  dynamics.  The  incorporation  of 
n  x  n  feedforward  transfer  function  L(s)  is  included  to  permit  compensation  for  load  distur¬ 
bances  which  can  be  sensed  prior  to  their  effect  being  observed  on  the  plant.  In  cases  where 
plant  delays  are  appreciable  this  feature  may  add  substantially  to  the  ability  to  compensate 
for  output  loads.  The  design  model  proposed  offers  several  significant  alternatives  for  con¬ 
trol  of  flexible  space  structure  control  where  a  variety  of  sensors  may  be  available  which  can 
predict  dynamic  loading  under  varying  environmental  conditions  such  as  thermal  gradients 
encountered  with  exposure  to  sun  light. 

Exogenous  disturbances  d,  sensor  noises  and  desired  set  points  u  are  modeled  as 

zero  mean,  wide  sense  stationary  random  processes  and  are  therefore  characterized  by  their 
respective  power  spectral  densities;  Gd(s),  Gn^s),  Gn/(s),  Gu(s).  Notice  that  no  assumption 
is  made  about  the  form  of  the  controller  C(s)  except  the  obvious  dimensions:  m  x  n. 

A  significant  feature  of  the  Youla  method  is  the  algebraic  parametrization  of  all  stabi¬ 
lizing  compensators  for  a  given  plant  and  sensor  combination.  Following  [3]  we  will  call  the 
combination  of  P  and  F  admissible  [3]  if  each  transfer  function  is  free  of  hidden  modes  and 
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Figure  2.1:  Multiloop  Control  Configuration  for  Wiener- Hopf  Design 

if  the  combination  does  not  introduce  cancellation  of  poles  and  zeros  in  the  Closed  Right 
Half  Plane  (C+).  (This  amounts  to  a  requirement  for  stabilizability  and  detectability  of  the 
plant  model.)  Given  that  P  and  F  are  admissible  let 

F(s)P(s)  =  D;'(s)N,(s)  =  N,(*)D;'  m 

where  Ni,  Dc  (resp.  Nr,  Dr)  are  left  (right)  coprime  factorizations.  As  a  result  there  will 
always  exist  a  pair  of  polynomial  matrices  X($)  and  Y’(,s)  such  that  [12,  4] 

Dc{s)X{s)  +  N6{s)Y{s)  =  L 

Then  the  closed  loop  system  is  asymptotically  stable  if  and  only  if  the  controller  is  of  the 
form 

C(s)  =  [Y(s)  +  Dr(s)K(s)}[X(s)  -  Nr(s)K(s)}~1 .  (2.41) 

where  K  is  any  m  x  n  real  rational  matrix,  analytic  in  C+  and  det[X(^)  —  Nr(s)K(s)}  0. 
Thus  all  stabilizing  controllers  can  be  expressed  in  terms  of  a  “stable  parameter11,  K{&). 

With  the  above  parameterization  in  hand  synthesis  can  be  based  on  an  optimal  con¬ 
trol  problem  which  attempts  to  minimize  tracking  error  e  =  u  —  y  subject  to  a  power-like 
constraint  on  the  control  r.  Thus  let1 

Jt  —  ~ — :E  jjY  J  e(^)e*(5)d^|  (2.42) 

subject  to  a  constraint 

J,  =  {Tr  J  PJ(3)r(j)r.(3)P,.(a)d3  j  ,  (2.43) 


xWe  use  the  notation  u*(s)  =  ur(~5),  E  denotes  expectation,  and  Tr  denotes  trace. 
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where  P,(,s)  is  the  transfer  function  from  the  control  signal  r  tp  those  sensitive  plant  modes 
which  must  be  protected  against  saturation.  The  resulting  optimal  control  problem  is  to 
find  the  controller  (i.e.  the  parameter  K(s))  such  that  the  cost 


J  —  Jt  +  kj9  (2.44) 

is  minimized  where  K(,s)  is  allowed  to  vary  over  all  stable  transfer  functions  of  appropriate 
dimensions  and  k  >  0  is  a  real  parameter  which  plays  the  role  of  a  Lagrange  multiplier  and 
effectively  permits  tradeoff  between  tracking  and  saturation. 

A  complete  solution  can  be  obtained  by  the  following  procedure: 

Wiener-Hopf  Design  Procedure  (Youla’s  Method) 

Step  1  Obtain  coprime  factorizations 

F(,)P(s)  =  D;'(s)N,(s)  =  K(s)d;'(s) 

Step  2  Compute  spectral  factorizations 

Dr.[P.P  +  kQ}Dr  =  A.A 
DtGeDt.  =  Q  Q. 

where  fi,  A  are  analytic  together  with  their  inverses  in  CRHP  and 

Ge  =  Gu  +  ( FP0  +  L)Gd{FP0  +  L).  +  F0Gn,F0.  +  L0GntL0., 

is  the  combined  power  spectral  density  of  the  exogenous  inputs  to  the  control 
2 

Step  3  Find  a  solution  X(s)F(s)  to  the  Diophantine  equation; 

zMa)X(s)  +  =  /p. 

Step  4  Compute  the  transfer  function 

Z(s)  =  Dr.P.[Gu  +  P0Gd{FPo  +  L).]D(, . 

and  the  stable  parameter2 3 

K(s)  =  a-1  ({A:1zfi:1}+  +  {ad;1}^},)  fr1  -  d^y.  (2.51) 

Then  the  optimal  controller  has  transfer  function  given  by  the  parametric  formula 

(2.41). 

2We  remark  that  the  computations  of  the  spectral  factors  in  (2.46)  (resp.  (2.47))  effectively  replace 
the  computational  step  of  solving  a  Riccati  matrix  equatiou  for  the  control  (resp.  filter)  problem  typically 
encountered  in  modern  state-space  methods  for  control  design.  For  systems  with  distributed  parameter 
models  the  Riccati  equation  is  a  PDE  and  therefore  difficult  to  solve  by  standard  methods. 

JA  rational  (matrix)  function  has  a  (partial  fraction)  expansion  A(s)  =  {A(s)}+  +  {i4(s)}_  +  {i4(s)}^ 
where  {.}+  (resp.  {.}-)  is  the  part  analytic  in  9?es  >  0-the  causal  part  (3Res  <  0-tlie  anti-causal  part)  and 
{.}x  is  the  part  associated  with  poles  at  infinity. 


(2.49) 

(2.50) 


(2.46) 

(2.47) 


(2.48) 

system. 


(2.45) 
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The  design  method  permits  direct  evaluation  of  control  design  tradeoffs  in  terms  of  certain 
critical  frequency  response  functions.  From  the  control  architecture  shown  in  Fig.  2.1  we  can 
readily  derive  the  closed  loop  frequency  responses  as 


y  =  PCS(u  —  F0nj  —  L0ni)  +  [ P0  —  (F  P0  +  L)CS\d  (2.52) 

r  =  CS[u  -  F0nf  -  L0nc  —  (FP0  +  L)d]  (2.53) 

e  =  (I  -  PCS)u  +  PCS(F0nf  +  L0nt)-[Po-  PCS{FP0  +  L))d  (2.54) 

where  the  system  sensitivity  operator 

S  =  [I  +  FPC}-\  (2.55) 


is  a  critical  transfer  function  whose  frequency  response  highlights  the  broadband  action  of 
the  feedback  control  law.  For  the  optimal  control  computed  by  the  method  outlined  above 
the  system  sensitivity  is  given  by 


S  =  {X  -  NrK)Dt , 


(2.56) 


a  relation  whose  simplicity  underlies  the  algebraic  characterization  of  feedback  design  from 
the  frequency  response  viewpoint. 

The  basis  for  resolving  the  engineering  tradeoff  of  closed  loop  control  action  in  terms 
of  the  conflicting  requirements  for  stability,  disturbance  rejection,  command  tracking,  and 
saturation  avoidance  is  the  application  of  Wiener-Hopf  method  for  the  solution  of  the  op¬ 
timization  problem  outlined  above.  It  can  be  easily  shown  that  the  optimal  performance 
objectives  Jt  and  J,  are  readily  computed  from  the  frequency  domain  solution  above  from 
the  following  formulae: 


where 


Ga 

Pd 

R 


Jt  = 

1  rj°° 

—  /  Tr  $t(ju>)du>, 

^7TJ  J  —joo 

(2.57) 

1  ri°° 

J ,  = 

—  Tr  $,(ju)du, 

Z'K ]  J -joo 

(2.58) 

(/  -  PR)GU(I  -  PR).  +  (. PR)Gi(PR ).. 

(2.59) 

QRG'R ., 

(2.60) 

FaGnf  F„.  -f  L0GneL0 

(2.61) 

FPo  +  L, 

(2.62) 

(Y  +  DrK)D(. 

(2.63) 

Thus  the  design  method  proceeds  by  resolving  the  tradeoff  between  tracking  performance 
(requiring  high  loop  gains)  and  saturation  limitation  of  the  critical  plant  modes  (which  will 
limit  loop  gains)  by  direct  implementation  of  the  composite  cost  relative  scale  factor  k  in 
(2.44)  which  plays  the  role  of  a  Lagrange  multiplier. 
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Remarks  on  the  Extension  of  Youla’s  Methods.  Recent  activity  in  the  theory  of 
robust  control  system  design  methods  has  focused  on  the  solution  of  certain  minimax  or 
worst,  case  design  problems  described  in  terms  of  frequency  response  models  and  frequency 
dependent  model  uncertainty  [12].  Solutions  have  been  obtained  for  a  class  of  H°°  optimal 
control  design  problems  by  extending  the  application  of  the  algebraic  characterization  of 
the  class  of  stabilizing  controllers  for  a  given  plant  given  in  (2.41).  In  Youla’s  original  work 
the  notion  of  coprime  factorization  was  based  on  the  natural  choice  for  rational  transfer 
functions:  the  factorization  was  with  respect  to  the  ring  of  polynomials.  This  paradigm 
ultimately  limited  the  popularity  of  Youla’s  original  design  methods  since  for  all  but  rather 
simple  systems  the  factorization  steps  became  computationally  cumbersome. 

The  modern  approach  described  by  Vidyassagar  [12]  extends  Youla's  characterization 
of  all  stabilizing  controllers  by  introducing  the  idea  of  obtaining  coprime  factorizations  of 
transfer  functions  over  the  ring  of  stable,  rational  transfer  functions.  It  is  shown  in  [12]  that 
the  set  of  all  transfer  functions  with  poles  contained  in  a  half  plane  form  an  algebraic  ring 
with  respect  to  the  usual  operations  of  addition  and  multiplication.  This  characterization 
is  of  practical  importance  for  control  system  design  since  one  can  now  further  constrain  the 
optimal  control  computation  so  that  a  prescribed  degree  of  stability  will  be  obtained  in  the 
closed  loop.  In  general,  the  extension  of  this  algebraic  characterization  of  the  class  of  sta¬ 
bilizing  controllers  for  irrational  transfer  functions  is  incomplete,  but  for  the  special  class  of 
irrational  transfer  functions  arising  in  flexible  structure  control  the  essential  characterization 
is  now  complete. 

2.3  Wiener-Hopf  Design  for  A  Class  of  Irrational  Transfer  Functions 

The  formal  extension  of  the  algebraic  computations  for  stabilizing  and  optimal  control  syn¬ 
thesis  in  the  frequency  domain  is  an  area  of  active  research  [12,  pp.  357].  In  this  section 
we  summarize  available  results  for  a  limited  class  of  ’pseudo-meromorphic*  transfer  func¬ 
tions  which  include  most  linear,  time-invariant  models  arising  in  control  of  flexible  space 
structures  [2].  We  refer  the  reader  to  recent  work  of  Baras  for  proofs  and  details  [13,  14]. 

The  results  of  Baras  extend  the  critical  constructions  required  for  Wiener-Hopf  design 
to  the  class  of  tpseudo-meromorphic,  transfer  functions,  denoted  as  containing  transfer 
functions 

GeH°°(  c+) 

where  C+  denotes  the  closed  right  half  plane  and  H°°  is  the  Hardy  space  of  bounded  analytic 
functions  on  C+.  Such  transfer  functions  have  (weak)  coprime  factorizations  over  M{A™)4 

T{iw)  =  Nr(iu;)D;l(iw)  =  D;l(-iu>)Nt{iu>) 

in  the  sense  that  Nr,  Dr ,  Dt ,  Nt  are  in  H°°( C_|_)  provided  that  Dt ,  Dr  are  inner  [4,  pp.  635]. 

Alternately,  we  say  that  N(s)(p  x  m),  D(s)(m  x  m)  are  strongly  coprime  if  there  exist 
8  >  0  such  that 

||Af($)x||  +  ||Z?(j)x||  >  8  >  0 

4  Let  M(A^i)  be  the  set  of  matrix  transfer  functions  (of  dimensions  determined  by  the  context)  with 
elements  in  A2. 
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for  all  x  G  Cn  and  s  G  C  +  .  Then  the  Carleson  Corona  theorem  asserts  that  N,  D  are  strongly 
coprime  if  and  only  if  there  exist  A",  Y  in  H™Xrn  such  that 

NY  +  DX  =  /m. 

Finally,  the  algebraic  constructions  required  for  the  Wiener-Hopf  control  design  method 
described  above  are  extended  by  the  following  result  [15,  16]. 

Theorem  1  Assume  that  P,  F,  L,C  G  A™  and  FP  G  A Then  all  of  the  following  hold. 

1.  The  singularities  of  [/  +  CFP\~lFP  are  the  zeros  of  det[D{X  +  NtY]. 

2.  There  exist  stabilizing  controllers  provided  that  Nt,Dt  are  strong  coprime . 

3.  If  Nr,  Dr  are  strongly  coprime  then  the  closed  loop  system  is  BIBO  (or  exponentially) 
stable  if  and  only  if 

C  =  {Y  +  DrK)(X-NrK)-1 

where  K  is  any  element  of  A™,  such  that  det(AT  —  NrK)  ^  0,  analytic  in  >  0 
(5?es  >  0),  and  X  and  Y  are  solutions  of  the  Diophantine  identities 

DtX  +  NtY  =  /. 

Computational  Requirements  for  Wiener-Hopf  Methods  For  rational  transfer  func¬ 
tions  the  required  computational  steps  outlined  by  Youla  and  his  colleagues  for  Wiener-Hopf 
control  can  all  be  carried  out  over  the  ring  of  polynomials.  More  recently,  the  approach 
outlined  by  Vidyassagar  [12]  suggests  that  by  performing  computations  over  a  ring  of  stable 
transfer  functions  (denoted  S  and  including  all  functions  rational,  proper,  and  analytic  in 
C+  U  {oo})  the  required  computations  can  be  directly  supported  in  terms  of  state  space 
realizations  for  the  individual  transfer  functions.  The  key  observation  is  that  the  collection 
of  all  proper,  rational  transfer  functions  with  poles  in  C+  form  an  algebraic  ring  which  has 
the  structure  of  a  principal  ideal  domain  (<S). 

Next  we  briefly  outline  several  options  for  extending  the  required  computations  to  ^4^. 
Just  as  in  the  rational  transfer  function  case,  the  nonuniqueness  of  state  space  realizations 
for  transfer  functions  may  offer  computational  alternatives.  Computational  requirements  for 
Wiener-Hopf  synthesis  include  the  following. 

1.  Coprime  Factorization.  For  rational  transfer  functions  the  coprimeness  condition  for 
a  pair  ( D )  essentially  amounts  to  a  requirement  for  noncollocation  of  zeros  for  N  and 
D.  In  general,  irrational  transfer  functions  may  not  have  coprime  factorizations [12]. 
In  our  studies  the  class  of  pseudo-meromorphic  transfer  functions  includes  most  mod¬ 
els  for  structural  vibration  control.  Thus  we  have  the  existence  of  (strong)  copriine 
factorizations  in  A™. 

2.  Solution  of  Diophantine  Relations.  Restricting  attention  to  pseudo-meromorphic 
transfer  functions  gives  the  required  existence  of  solutions  to  (2.49).  Despite  the  poten¬ 
tial  lack  of  convergence  of  Euclidean  division  in  *4^Berenstein  and  Struppa  [17]  have 
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obtained  explicit  solutions  to  equations  of  the  form  (2.49)  for  ( TV,  D)  strongly  copriine. 
The  formulae  obtained  in  [17]  are  nonalgebraic  and  extremely  complex.  Approximate 
solutions  to  the  Bezout  relations  can  also  be  obtained  by  sampling  the  Fourier  response 
of  the  individual  terms  Z?*,  Nt>X,Y  on  a  finite  set  of  frequency  samples.  This  reduces 
the  computational  problem  to  that  of  solving  a  set  of  symultaneous  linear  equations. 

The  insight  from  the  algebraic  approach  to  synthesis  of  stabilizing  control  [12]  indicates 
that  a  particular  solution  X,Y  to  (2.49)  is  a  nominal  stabilizing  feedback  control ;  K  = 
YX~l.  For  most  models  of  flexible  structures  arising  in  space  applications  the  plant 
is  nominally  stable  except  for  a  finite  number  of  poles  either  at  the  origin  (associated 
with  rigid  body  inertias)  or  in  the  right  half  plane.  For  such  systems  we  have  obtained 
direct  methods  for  computing  the  required  (stable)  coprime  factorizations.  (See  for 
example  benchmark  problem  2.) 

3.  Spectral  Factorization  and  Causal  Projection.  In  most  cases  arising  in  flexible 
structure  control  the  resulting  irrational  transfer  function  model  is  meromorphic.  In 
some  simple  cases  the  computational  problems  can  often  be  reduced  by  the  identifica¬ 
tion  of  an  alternate  state  realization  [18].  A  more  general  approach  is  to  use  frequency 
sampling.  Methods  for  matrix  spectral  factorization  have  been  studied  as  alternatives 
to  solving  infinite  dimensional  Riccati  equations  by  Davis  [19]  and  Bennett  [20]  where 
a  numerical  algorithm  is  described.  In  the  next  section  we  describe  results  of  testing  a 
spectral  factorization  algorithm  based  on  frequency  response  sampling. 

4.  Linear  Fractional  Combinations.  All  remaining  constructions  involve  linear  frac¬ 
tional  transformations  which  can  be  readily  supported  via  either  explicit  transfer  func¬ 
tion  computations  or  frequency  sampling. 

Our  conclusion  from  this  survey  is  that  frequency  sampling  offers  an  attractive  alternative 
for  the  computational  requirements  of  Wiener-Hopf  synthesis  and  can  readily  be  extended  to 
transfer  functions  of  pseudo- meromorphic  type.  The  computational  approach  then  obtains  a 
frequency  sampled  approximation  to  the  ideal,  optimal  controller  for  the  given  design  problem 
and  essentially  provides  a  specification  for  such  control.  We  remark  that  in  contrast  to 
state-space  computations  this  approach  decouples  the  realization  of  the  controller  from  the 
computation  of  its  response. 

3  Phase  1  Project  Results:  Computation  of  optimal  control  for 
flexible  structures  by  frequency  sampling. 

3.1  Frequency  Response  Design  Method  for  Vibration  Control 

In  the  background  material  of  section  2  we  have  described  the  basis  for  the  design  approach 
we  have  investigated  in  Phase  1  study  for  flexible  space  structures.  Wiener-Hopf  optimiza¬ 
tion  is  the  basis  for  evaluating  engineering  tradeoffs  with  respect  to  the  choice  of  required 
control  law.  For  LQG  type  designs  the  tradeoff  analysis  from  quadratic  optimization  results 
in  a  pair  of  feedback  gain  matrices.  In  the  method  we  employ,  the  result  is  a  complete 
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controller  frequency  response  model.  The  Phase  1  effort  has  focused  on  the  demonstration 
of  the  required  frequency  response  modeling  and  control  design  computations  for  some  sim¬ 
ple  benchmark  problems  in  vibration  control  and  isolation  for  flexible  structures.  We  have 
demonstrated  the  feasibility  of  solving  the  Wiener-Hopf  problem  using  frequency  sampled 
data  representation  for: 

1.  the  plant  model  transfer  function,  including  the  frequency  response  of  the  elastic,  struc¬ 
ture  to  excitation  by  localized  (point  wise)  actuation  and  possibly  distributed  distur¬ 
bance  excitation, 

2.  the  sensor/actuator  model  transfer  function,  including  provisions  for  both  feedforward 
and  feedback  processing, 

3.  the  PSD  of  exogenous  driving  signals,  including  sensor  noises,  output  disturbances,  and 
command  inputs, 

4.  control  objectives,  including  the  specification  of  desired  closed  loop  performance  in 
terms  of  frequency  dependent  quadratic  objective  functions  for  both  tracking  and  sat¬ 
uration  avoidance. 

In  our  benchmark  problems  we  have  considered  both  irrational  transfer  functions  arising 
from  standard  model  of  a  flexible  beam  and  finite  element  model  of  a  hybrid  structure 
representing  a  problem  of  active  vibration  isolation.  The  method  performs  with  success  in 
both  cases.  Central  to  the  success  of  the  effort  is  a  set  of  computational  algorithms  and  an 
interactive  environment  for  computation  with  frequency  response  models.  One  feature  of 
the  sampled  frequency  response  computations  is  the  requirement  to  choose  a  bandwidth  and 
sampling  interval  for  the  model  representation.  Although  guidance  for  these  choices  can  be 
obtained  from  analysis  of  the  transfer  function  models,  it  has  been  our  experience  that  some 
fine  tuning  is  desirable  in  applications.  The  interactive  environment  for  frequency  sampled 
computations  supports  this  requirement. 

3.2  Computational  Algorithms  for  Frequency  Domain  Control  Design 
3.2.1  Coprime  (stable)  Factorization 

Many  flexible  structure  models  involving  hybrid  dynamic  models  often  include  a  finite  num¬ 
ber  of  unstable  poles.  In  such  cases  the  coprime  factorization  of  the  plant  transfer  function 
is  a  critical  step  which  will  guarantee  the  stability  of  the  closed  loop  control  system.  The 
importance  of  the  class  of  pseudo-meromorphic  transfer  functions  is  that  the  notion  of  stable 
copnme  factorization  can  be  extended  to  these  systems. 

By  way  of  review,  recall  that  a  rational  p  x  m  strict  proper  transfer  function,  P(s)  has  a 
minimal  realization  in  terms  of  state  space,  matrix  parameters  (A,  B,  C) 

P(s)  =  C(sl  -  A)~lB,  (3.1) 

where  ( A,  B )  is  controllable  and  (A,C)  is  observable.  Then  the  notion  of  coprime  factor¬ 
ization  with  respect  to  a  class  of  “stable”  systems  for  which  all  members  have  poles  within 
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an  arbitrary  (left)  half  plane,  can  be  connected  with  state  space  computations  as  follows. 
Let  5  be  the  class  of  transfer  functions  with  poles  p,-  such  that  We  p,-  <  a.  Given  a  re¬ 
alization  (A,B,C)  for  P(s)  then  under  the  above  conditions  there  exist  5-stable  coprime 
factorizations, 

P  =  DJ' Nt  =  NrD;\  (3.2) 

which  satisfy  a  doubly  Diophantine  relation, 


* 
_ 1 

'  Dr  ~Xt  ' 

1 

s 

b 

* 

(3-3) 


It  can  be  easily  verifed  that  the  individual  transfer  functions  Nr,  Nt,  Dr,  D(,  Xr ,  Yr,  N(,  D(  G 
5  can  be  obtained  as: 


Nt  =  C(sl  —  A„)~XB, 

Dt  =  I  —  C(sl  -  A0)~1F, 

Nr  =  C{sl  -  Ac)~lB, 

Dr  =  I  -  H(sl  -  AC)~1B,  (3.4) 

*r  =  HisI-AJ-'F, 

Yr  =  1+  H{sl  -  A0)-'B, 

Xt  =  H(sl  —  AC)~1F, 

Yt  =  I  +  C(sI-Ac)~lF , 

where  the  n  x  n  matrices  Ac  —  A- f  BF  and  A0  =  A  - f  HC  are  constructed  by  appropriate 
choice  of  the  state  feedback  gain  matrix  F  and  the  output  injection  matrix  H  so  that  the 
eigenvalues  of  Ac  and  A0  are  contained  in  S. 

The  computations  of  coprime  factorization  can  be  readily  extended  to  the  class  of  models 
representing  the  dynamics  of  flexible  structures  under  the  following  assumption.  Let  P(s) 
be  an  irrational  transfer  function  such  that 

P(,)  =  Ps(s)  +  Pi(»)  (3.5) 

where  Ps  G  5  C  and  Pj  is  rational  and  has  (a  finite  number)  of)  poles  outside  the  half 
plane  defining  5;  i.e.,  P(^)  has  only  a  finite  number  of  “unstable”  poles.  In  this  case  we 
can  obtain  the  stable  coprime  factorization  for  P(s)  with  respect  to  5  as  follows.  Obtain  a 
coprime  factorization 

ps  =  NrD;\ 

the  unstable,  rational  term.  This  can  be  achieved  by  the  state  space  constructions  described 
above  without  reference  to  Ps — the  irrational  spectra  of  P.  Then  P  has  coprime  factorization 

P  =  NrD 71  =  [Nr  -  PsDr)b;\  (3.6) 

where  Nr,  Dr  are  5-stable.  The  separation  of  terms  in  (3.5)  is  readily  carried  out  given  P(s) 
by  computing  the  residues  of  the  finite  number  of  unstable  poles  contributing  to  Pj.  For  most 
flexible  structure  problems  the  number  of  unstable  poles  is  small  arising  from  the  interaction 
of  elastic  dynamics  of  the  distributed  structure  with  localized  rigid  body  dynamics. 
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3.2.2  Solution  of  Diophantine  Relations 

The  characterization  of  the  class  of  all  stabilizing  controllers  for  a  given  rational  transfer 
function,  P(.s),  as  obtained  by  Youla  and  his  coworkers,  is  given  in  terms  of  a  coprime 
factorization  P  =  D^1  and  solution  to  a  Diophantine  relation  (2.49).  All  stabilizing 
controllers  C(s)  for  P(.9)  can  then  be  described  by  the  linear  fractional  transformation 

C=  [X  -  KNt]~l[Y  -  KDt], 


The  characterization  is  algebraic  for  rational  transfer  functions  and  the  computations  can 
be  effectively  carried  out  by  state  space  computations  as  descirbed  in  the  previous  section. 

For  irrational  transfer  functions  for  which  (strong)  coprime  factorizations  exist  we  can 
readily  obtain  solutions  to  the  Diophantine  relation  by  frequency  sampling.  Rewriting  (2.49) 
in  the  form 


(£><(.),  AT, (,)) 


Y(>) 


=  4. 


(3.7) 


it.  is  clear  that  the  required  computation  obtains  a  right  inverse  for  the  p  x  m  +  p  matrix 
[Dt(s)}  N*(,s)].  This  is  readily  obtained  by  frequency  sampling  using  one  of  several  standard 
numerical  algorithms.  A  numerically  stable  approach  is  to  use  Singular  Value  Decomposition 
(SVD)  of  the  p  x  m  +  p  matrix  obtained  at  each  complex  frequency  sample  s.  The  routine 
CSVDC  contained  in  the  Linpack  software  provides  a  direct  way  to  obtain  the  right  inverse 
(i.e.,  Moore-Penrose  pseudo-inverse)  from  the  SVD. 


3.2.3  Spectral  Factorization  and  Causal  Projection 

A  critical  computational  requirement  of  the  design  method  we  have  investigated  is  the  re¬ 
quirement  for  solving  two  instances  of  Wiener-Hopf  optimization  based  on  frequency  response 
data.  Spectral  factorization  arises  in  optimal  control  and  filtering  problems  by  association 
with  the  solution  of  a  Wiener-Hopf  integral  equation; 

[  h(t  —  r)w(r)dr  =  f(t)  (3.8) 

Jo 

for  t  >  0  where,  for  example,  h(t)  may  be  the  covariance  matrix  of  some  noise  process,  f(t) 
is  a  specified  causal  function,  and  u>(£),  for  t  >  0  is  the  solution  sought.  The  Wiener-Hopf 
technique  solves  (3.8)  by  the  identification  of  a  factorization  of  the  Laplace  transform 

H(s)  =  r  h(t)e-3tdty 

J —00 

of  the  form 

H(s)  =  F(s)Ft(-s)  (3.9) 

with  the  property  that  F(s)  is  analytic  together  with  its  inverse  in  the  closed  right  half  plane, 
C+.  For  the  case  of  H(s)  a  rational  (possibly  matrix  valued)  transfer  function,  existence 
and  uniqueness  of  the  spectral  factorization  above  follows  under  the  conditions: 


1.  H(s)  —  H(s);  i.e.,  H(s)  is  the  transform  of  a  real- valued  function  h(t). 
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2.  H(s)  =  Hm(s)\  i.e.,  H(s)  is  “para-hermittian” 

3.  H(s)  is  of  normal  rank;  i.e.,  full  rank  almost  everywhere  in  C. 

4.  H(itu)  is  positive,  semi-definite  for  uj  €  K. 

More  generally,  (following  the  classical  theory  of  Gohberg  and  Krein  [21])  we  have  for  the 
typical  optimal  filtering  or  control  problem  that  H{s)  =  I  +  G«(s)G(s)  where  G(s)  is  the 
transfer  function  of  the  system  to  be  controlled  and  under  the  assumption  that  G*G  is: 

1.  positive,  semi-definite  for  s  =  iu, 

2.  G+G  is  the  transform  of  a  function  which  is  both  L\  and 
then  the  spectral  factorization  (3.9)  exists  and  has  the  property 

F(iw)-Ie  T{LtnL$)  (3.10) 

(where  L +  denotes  those  L\  functions  with  positive  support)  and  F(s)  =  F(i). 

We  are  interested  in  obtaining  a  consistent,  computationally  efficient  approximation  for 
spectral  factorization  problems  arising  from  transfer  functions  G(s)  belonging  to  the  Hardy 
space  H 2  fl  H°°  (which  may  include  irrational  cases).  The  approach  we  have  in  mind  involves 
sampling  and  interpolation  of  the  frequency  response  data  H(iuj ). 

We  will  also  consider  the  associated  problem  of  causal  projection ; 

V+  {/  +  r  f(t)e~ivtdt}  =  /  +  f°  f(t)e~iutdt.  (3.11) 

J-oo  JO 

Then,  in  the  scalar  case,  if  F(s)  is  the  scalar  caused  factor  in  (3.9)  and  $(a)  =  In  H(s)  we 
can  write 


*(»)  =  «»)}+  ■P-  {*(•)}  (3.12) 

=  In  F(.)  +  In  F.(s)  (3.13) 

so  that  the  spectral  factor  can  be  obtained  as  F(s)  =  exp  V+  {In  i/(s)}.  In  the  matrix  case  we 
employ  a  recursive  algorithm  for  spectral  factorization  which  is  loosely  based  on  a  Newton- 
Raphson  iteration  for  an  associated  Riccati  equation.  First,  we  consider  computation  of 
causal  projection  and  spectral  factorization  in  the  scalar  case  by  frequency  sampling. 

Spectral  Factorization  of  Scalar  Data  Sequences.  In  the  Phase  1  effort  we  have  devel¬ 
oped  and  tested  computer  algorithm  for  causal  spectral  factorization  of  frequency  response 
samples  based  on  an  interpolation  method  of  F.  Stenger  [22].  Roughly,  Stenger’s  method 
offers  an  approximation  of  a  function  g(w)  €  T  ( Lp )  a s  an  expansion 

5(a)(^)  =  £  9((j  + 

j--oo 


(3.14) 


SEI-89-03-15-WB 


20 


where  Au;  is  a  fixed  sampling  interval,  and  Xj  (<*’)  is  the  jth  (approximate)  characteristic 
function  which  approachs  the  ideal: 

x;M  = 

We  summarize  the  theoretical  basis  for  Stenger’s  method  in  the  following  theorems. 
Theorem  2  Let  g  G  T  ( Lv ).  Then 


1,  for  u>  G  +  l)Au>) 

0,  else 


(3.15) 


0(a)H=  £  9({j  +  ^)Au>)xj(u;) 


(3.16) 


approximates  g  in  the  sense  that  \\g  —  <^a)||p  — *  0  as  (Au/,  k)  — ►  (0+,  1  —  )  where  the  approxi¬ 
mate  characteristic  functions  are 


X;(“>)  =  -t=[!  +k<j>j{u)% 


and 


4>j(  w)  =  sn 


2  K 

—  log 

7 n 


2  y/k 

'u>  —  ( j  +  l)Au/ 


-K  +  iK\  k 


(3.17) 


(3.18) 


uj  —  JAoj 

Here  we  use  the  standard  notation  for  the  elliptic  functions  with  parameter  fc,  0  <  k  <  1; 


vizM 


if  and  only  if 


z  =  sn[u,  k] 


,  s  r  dt 

u{z)  =  /  / 

Jo  i/a  -t2)d  - 


v/(i  -  f2)(i  - fc2f2)’ 

A'(fc)  =  Ti(l),  and  A"  =  A'(Vl  -  fc2). 


(3.19) 

(3.20) 

We  remark  that  the  characteristic  function  Xji^)  given  in  (3.17)  has  the  following  prop¬ 
erties  [22]: 

1.  for  a;  €  (j Aw,  (j  +  l)Aw)  =►  X»  €  p1/2  +  fc"1/2)/2,  Hr*/*], 

2.  for  u;  G  9?  —  [j  Au;, (j  +  l)Au;]  =>  XjP)  €  [0,(A;-1/2  —  fc1/2)/ 2]  . 

Theorem  3  The  characteristic  function  Xjp)  has  the  explicit  representation 


m  =  —  oo  lu >  ~  Pm  UJ  -  pm  ) 

(3.21) 

where  the  poles 

Pm  =  ( J  +  am)Aui 

(3.22) 

and  residues 

™  m2 

"  "  iy/kK^  “ 

(3.23) 
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depend  on  the  parameters 


— 


1  —  iq 

«  =  e*p(-»fl  =  «P 


dt 


^(1 -«’)(! -(!-*’)«’) 


]• 


(3.24) 

(3.25) 


We  remark  that  (p^\rm)  corresponds  to  the  anti-causal  (right  half  plane)  poles  and 
residues  of  Xj  while  (p^\fm)  are  causal.  This  motivates  the  following  result. 

Theorem  4  Given  g(u>)  €  T  ( Lp )  the  causal  projection  of  g(w),  V+  {<7},  can  be  obtained 
in  the  limit  as  follows, 


9+]  =  / *  ..lin,1nx  ,  ,  Y  3{{j  +  {.XjM} 

(Au>,fc)-.(0+,l-)  j__00  2. 


lim 


Y  9(U  +  o)Aa;)  Y 


(Au,,feH(0+,i-)j.^oo''"'  2'  '  ^  ... 


(3.26) 


m  =  —  oo  UJ  Pm 

Then  ||^a)(u;)  -  V+  {<7(u;)}||p  -  0  as  ( Au/,  k )  -4  (0+,  1-). 

Now  assume  that  the  scalar  function  h( u>)  €  T  ( L\  (~l  )  is  positive  real  and  thus  has 

a  spectral  factorization  h(io)  =  f(us)f(uj)  where  /( u>)  is  the  unique  causal  factor.  We  can 
extend  the  above  causal  projection  computation  vja  logarithmic  transformation  of  the  data 
as  follows. 

Theorem  5  Let  h(s)  have  the  above  properties,  then  the  spectral  factorization  /i(s)  = 
f.(s)f(s)  exists  and  is  unique  where  f(s)  and  1  / f(s)  are  both  analytic  in  C+.  Then  the 
approximation 


/(o)(u>)  =  exp[  Y  lnM0'  +  oAu,)^+  {Xi(^)}| 


(3.27) 


J=-00 


with  V+  {Xj^)}  given  as  in  (3.26)  has  the  property  ||/^(u>)  —  /(u>)||p  — ►  0  for  both  p  =  1,2 
as  ( Au>,  k)  — ►  (0+,  1  —  )• 

By  way  of  summary  of  the  theoretical  basis  for  the  algorithm  we  have  in  mind  we  direct 
attention  to  the  following  formula  for  (approximate)  causal  projection: 


V+  *  Y  9((j  +  \)±»)  Y 

j=-oo 


m  =  —  oo  Pm 


(3.28) 


and  for  (approximate)  causal  spectral  factorization: 

/M  «  exp[  Y  lnh((j  +  |)Aa;)  Y  '"'"(j)]  (3-29) 

j—  —  oo  "  m  =  — oo  W  Pm 

where  the  sampling  interval  Au>  is  small  enough  and  k  <  1  but  close  to  1.  Substitution  of 
the  expressions  (3.22)-(3.25)  reveals  the  relation 

fm  -7 rqma2 


u>  —  4 y/kK  [(u>  —  (ji  +  |)Au>)  -f  |  +  am] 


(3.30) 
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so  that  an  Mth  order  approximation  to  the  projected  characteristic  function  for  sampling 
interval  Au>  is 

M  f  , 

=  E  — ^  =•  V*(»  -  (J  +  5)Au,)  (3.31) 

so  that  the  computation  of  (3.28)  (or  (3.29))  for  uniform  sampling  io  =  (Au>,  t  =  0,  ±1,  ±2, . .  ., 
reduces  to  the  discrete  convolution  of  the  data  sequences  <7((,7  +  ^)Au;)  (resp.  In  h((j+  ^)Au;)) 
and  y>+( uj  —  ( j  +  |)Au;)  for  j  =  0,  ±1,  ±2, . . ..  This  is  efficiently  implemented  using  FFT 
processing. 

The  summation  (3.31)  converges  rapidly  if  the  characteristic  function  is  chosen  carefully; 
i.e.,  if  k  <  1  and  close  to  1.  Well  known  properties  of  elliptic  integrals  suggest  that  we 
compute  (k,K)  by  first  choosing  q  <  1,  then 

*=(!)2  “d 

where 

«.  =  2  E 

n  =  l 

^3  —  1  +  2  qn  . 

n=l 

For  example,  the  choice  q  =  \  obtains  k  =  0.999994  and  K  =  7.11943.  In  this  case  we  found 
M  =  8  sufficient  for  single  precision  computation  on  a  VAX  11/750. 

Spectral  Factorization  of  Matrix  Frequency  Samples.  For  the  general  MIMO  design 
problem  typical  of  flexible  space  structure  control  problems  the  spectral  factorization  steps 
are  required  for  matrix  valued  transfer  functions.  In  Phase  1  effort  we  have  implemented  and 
tested  a  computer  algorithm  which  operates  on  the  frequency  sampled  data  by  a  recursive 
procedure  to  obtain  the  spectral  factor.  The  basis  for  the  algorithm  can  be  obtained  by 
association  of  the  spectral  factorization  with  the  solution  of  a  Riccati  equation  arising  in  the 
context  of  a  quadratic  optimal  control  problem: 

min/*  |M0||2  +  \\y(t)\\2dt  (3.34) 

U 6M,„|  Jo 

subject  to  the  linear,  time-invariant  system  model 

i(t)  =  Ax(t)  +  Bu(t),  i(0)  =  a*  (3.35) 

y(t)  =  Cx(t),  t  >  0,  (3.36) 

where  we  assume  [A,  B,  C ]  is  a  minimal  realization  for  the  transfer  function  C7(s)  =  C[sl  — 
A\~lB.  The  optimal  control  is  known  to  be  a  linear  state  feedback  u(t)  =  —Koptx(t)  = 
—  BTPx(t)  where  P  is  the  unique,  positive  definite  symmetric  solution  to  algebraic  (matrix) 
Riccati  equation, 


(3.32) 

(3.33) 


PA  +  AtP  -  PBtBP  +  CTC  =  0. 


(3.37) 


SEI-89-03-15-WB 


23 


Standard  algebraic  manipulations  based  on  (3.34 )— (3.37)  provide  the  spectral  factorization 
relation  [23,  pp.  68]  which  we  write 

H{s)  =  1  +  Gt(-s)G(s)  =  Ft{-s)F{s)  (3.38) 

where  F(s)  =  I  +  K^lsI  —  A]-1  B  is  the  causal  spectral  factor  of  the  positive  real  transfer 
function,  H(s). 

The  algorithm  we  employ  was  first  suggested  by  Davis  and  Dickinson  [19]  and  takes  the 
form  of  a  Newton-Raphson  recursion  for  the  spectral  factor; 

Fn+.M  :=  V+  {[Fn-(ia,)l-1ff(iu,)[F„(iu;)|-‘}  F„(iw),  (3.39) 

where  V+  is  the  causal  projection  operator.  The  recursion  (3.39)  can  be  implemented  in  a 
form  which  enhances  its  numerical  properties  and  provides  an  effective  computer  algorithm; 

[Fn+1]-'  :=  [FJ-1  (/  +  F+  {[F^]-1  H  [F.r1  -  /})"’ .  (3.40) 

By  initializing  with  Fo  (an  m  x  m  diagonal  matrix)  with  diagonal  elements  equal  to  the 
spectral  factors  of  the  diagonal  elements  of  H  the  second  term  of  (3.40)  remains  a  pertur¬ 
bation  of  the  identity  (since  [F^]-1/f[Fn]-1  —  /  — +  0)  which  regularizes  the  computations. 
The  diagonal  initialization  guarantees  that  the  first  residual  [Fo]-1/f[Fo]-1  has  ones  on  the 
diagonal  and  all  off  diagonal  elements  less  than  one  in  magnitude.  The  resulting  numerical 
problem  is  well-conditioned. 

The  conceptual  algorithm  can  be  summarized  in  the  following  pseudo-code. 

Algorithm  for  Matrix  Causal  Spectral  Factorization 

Given:  An  m  x  m  matrix  valued  data  sequence  H(w)  defined  on  a  set  of  uniformly  sampled 
frequency  points  u ;  6  fl  together  with  the  property  that  at  each  u>  the  matrix  H(u> )  is  positive 
semi-definite,  symmetric. 

Initialize: 

F(0)(u/)  :=  diag{/fc(u/),  k  =  1, . . . ,  m} 

where  the  scalar  data  sequences  /fc(u>)  are  computed  using  (3.29)  from  the  m  diagonal  data 
sequences  of  H(u). 

[F-1(u>)](0)  :=  diag{l//fc(u/),A:  =  1,. . .  ,m}. 

Repeat  Until  (  ||i?^(u>)||  <  e)  where  e  is  a  specifed  tolerance. 

j  :=  j  +  1 

[F-'HP  :=  [F- +  K 
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Issues  in  Design  of  Efficient  Computer  Code.  The  primary  observation  that  permits 
the  design  of  an  efficient  computer  algorithm  follows  from  the  relation  (3.31)  so  that,  in 
effect.,  (3.28)  or  (3.29)  are  discrete  convolution  of  infinite  data  sequences.  In  applications, 
the  summations  appearing  in  (3.28)  or  (3.29)  will  converge  rapidly  and  can  be  appproximated 
within  any  desired  precision  by  finite  sums.  For  example,  in  design  of  control  laws  for  flexible 
structures  it  is  assumed  that  the  transfer  function  model  for  the  structural  flexure  response 
to  control  actuation  is  strictly  proper  with  effective  system  bandwidth  u>bw  which  is  known 
a  priori  (even  if  the  exact  transfer  function  is  not  known).  Modeling  issues  relative  to  the 
determination  of  u>bw  for  control  of  flexible  structures  are  considered  in  [24].  Once  the  model 
bandwidth  ujbw  has  been  obtained,  together  with  a  choice  of  a  frequency  sampling  interval 
Au>  then  we  compute  the  approximating  sequence  V+  {^(w)}  as 

Np 

P+  =  £  9{{j  +  J)AwV+(w  -  (j  +  i)Aw),  (3-41) 

:=-Nr 

where  Np  =  <jl>bw •  This  is  a  convolution  of  two  data  sequences  g,  of  length  2 Np. 

Since  in  many  applications  these  data  sequences  may  be  relatively  long,  an  efficient 
implementation  of  the  required  convolution  may  be  obtained  using  an  FFT  algorithm.  The 
development  of  reliable  computer  code  for  FFT  processing  has  been  quite  extensive  and 
several  efficient,  portable,  and  reliable  codes  are  available  in  the  public  domain  [25].  We 
have  found  the  Fortran  routine  FFT842,  available  in  [25],  an  efficient  alternative  which 
offers  support  for  the  required  complex  valued,  frequency  response  data  sequences. 

In  some  cases,  we  may  wish  to  operate  directly  on  data  sequences  obtained  from  mea¬ 
surements  on  physical  systems.  For  such  cases  we  must  concern  ourselves  with  potential 
numerical  errors  due  to  aliasing  and  Gibbs  phenomenon  [26].  For  application  to  control 
system  design  we  believe  it  is  important  to  consider  such  applications  in  the  development 
of  a  comprehensive  computer  code  and  we  have  included  various  standard  data  windows  for 
weighting  the  data  sequences  for  FFT  based  convolution. 

One  critical  tradeoff  in  the  development  and  testing  of  computer  code  to  support  the  re¬ 
quired  computations  for  Wiener-Hopf  control  system  design  is  storage  requirements  for  the 
complex  data  sequences.  In  Phase  1  the  focus  has  been  on  the  development  of  a  prototype 
computer  code  for  testing  the  numerical  algorithms  as  well  as  experimenting  with  the  choice 
of  bandwidth  and  data  sampling.  Thus  we  have  made  provisions  for  graphical  plotting  of 
the  various  data  sequences  as  the  algorithm  procedes.  We  have  made  explicit  provisions  for 
storing  a  number  of  data  sequences  which  are  useful  primarily  for  debugging  and  testing 
the  code.  As  a  result  the  prototype  code  is  a  data  storage  intensive  implementation  of  the 
method.  One  area  where  data  storage  can  be  efficiently  used  is  in  the  implementation  of 
the  spectral  factorization  by  consideration  of  several  basic  properties  of  the  matrix  data  se¬ 
quences.  Since  by  assumption,  the  model  (3.35)— (3.36)  involves  real,  vector-valued  functions, 
«(£),  t/(£)>  and  z(£),  the  transfer  function  matrix  has  the  property  G(iu>)  =  G(  —  iu>).  By  the 
construction  (3.38)  H(u>)  is  Hermittian,  positive  definite  for  each  u. 

At  the  jth  iteration,  the  recursion  is  driven  by  a  matrix,  residual  data  sequence; 

RU){ o>)  :=  FfT{-w)H{u)Fr'(w)  -  I 


(3.42) 
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Matrix  Sequence 

X{u)  =  X{-w) 

Hermittian 

G 

X 

H 

X 

X 

R 

X 

X 

V+{R} 

X 

F 

X 

Table  3.1:  Properties  of  Matrix  Data  Sequences 


which  is  also  Hermittian  for  each  u ;.  Table  3.1  summarizes  the  properties  of  the  matrix  data 
sequences  encountered. 

The  fact  that  the  data  sequences  involve  transforms  of  real-valued  processes  permits 
computational  savings  of  the  order  of  2NP  for  each  scalar  sequence  since,  at  the  jth  iteration, 
we  can  obtain  the  required  sequences  by  only  computing  Fj"1(u;)  for  >  0.  Let 

R$( u>)  denote  the  k  x  f  entry  of  the  residual  sequence  at  iteration  j.  Then  the  following 
holds; 

f>+ {**?(«)  =  «i.<V) - -P-  {rtuW)} 

=  *i j>(»)  -  V*  {*g(-u>)}  (3.43) 

=  fl!iV)  -  n  {«{?(«)}, 

where  we  have  used  the  facts  that  for  any  x(u>)  with  x{uj)  =  x(— u>),  V L  {x(u>)}  =  P+  {x(—u;)}, 
and  by  construction  of  the  residual  matrix, 

Rkt{—w)  —  Rki(w)  =  Rik(w). 

Thus  we  can  construct  the  lower  triangular  part  of  directly  from  its  upper 

triangular  part.  If  R(J\u>)  is  N  x  N  then  we  are  required  to  compute  causal  projection  using 
the  convolution  (3.41)  for  only  scalar  data  sequences  to  construct  the  required  causal 

projection. 

Finally,  the  computation  (3.40)  requires  symultaneous  solution  of  a  set  of  N  linear  equa¬ 
tions  for  each  to.  This  can  be  implemented  efficiently  using  the  routines  CGESL  and  GGEOO 
available  in  Linpack  [27], 

We  have  coded  and  tested  the  spectral  factorization  algorithm  for  several  problems  of 
both  scalar  and  matrix  types  and  including  both  rational  and  irrational  spectrum.  In  the 
next  section  we  described  the  results  of  these  studies  in  detail.  It  is  informative  to  monitor 
the  convergence  of  the  recursion  (3.40)  by  examining  the  residual  at  each  iteration  in  terms 
of  two  norms: 

max_ele(i?^)  :=  max{i?|/?(c<;)} 

a )tktl 

max.tr :=  max{<r(^(/)(a;)il^)(a;))}. 
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START  UP  MENU  AND  STATUS 

CHOOSE  A  NUMBER  BELOW. 

CURRENT  VALUE 

1.  INPUT  BAND  WIDTH. 

100.00 

2.  INPUT  NUMBER  OF  SAMPLE  POINTS. 

3.  PLOT  P(S) . 

4.  START  PROGRAM. 

5.  RETURN  TO  MAIN  MENU. 

1024 

Figure  3.1:  An  example  of  the  start.up  menu. 

3.3  Development  and  testing  the  interactive  software  environment 
3.3.1  Objectives  of  the  prototype  software  environment 

There  were  four  main  objectives  to  the  design  of  the  software  environment:  to  support 
the  frequency  sampling  analysis,  to  monitor  the  causal  spectral  factorization  algorithm  con¬ 
vergence,  to  resolve  design  trade-offs,  and  to  provide  a  means  to  obtain  desired  output. 
Additionally,  an  efficient  user  interface  consisting  of  a  menu  format  with  options  grouped 
according  to  major  functions  was  desired. 

Frequency  sampling  analysis  is  concerned  with  evaluating  the  frequency  responses  of 
specific  functions  at  discreet  frequency  values.  The  software  environment  was  designed  to 
perform  uniform  frequency  sampling  computations  using  a  maximum  of  1024  points.  The 
program  permits  the  user  to  select  a  bandwidth  and  a  number  of  frequency  samples  for  the 
computations  in  the  start.up  menu.  In  addition,  after  making  these  selections,  the  user  can 
plot  the  plant  transfer  function,  P(s),  to  evaluate  these  choices.  The  start.up  menu  appears 
in  figure  3.1. 

Another  objective  of  the  software  design  was  to  allow  the  user  to  monitor  the  convergence 
of  the  spectral  factorization  algorithm.  This  algorithm  uses  an  iterative  approach  to  converge 
on  the  spectral  factor.  The  results  of  the  convergence,  which  are  the  maximum  trace  and 
maximum  element  of  the  residual  and  the  number  of  iterations  needed  to  converge,  are 
dislpayed  in  a  summary  format,  as  seen  in  figure  3.2.  The  user  can  change  the  convergence 
threshold,  epsilon,  and  the  maximum  number  of  iterations  the  program  will  perform  before 
reaching  epsilon,  and  then  recompute  the  spectral  factor  if  desired. 

An  important  aspect  in  the  design  approach  of  this  paper  is  the  design  tradeoff  between 
tracking  cost  and  saturation  cost.  The  software  environment  enables  the  user  to  generate  a 
curve  of  saturation  cost  versus  tracking  cost  as  a  function  of  the  parameter  k.  This  curve 
can  be  seen  in  figure  3.3.  With  this  plot,  the  user  can  choose  a  k  value  that  corresponds  to 
some  maximum  saturation  level  and  a  corresponding  tracking  cost,  and  then  generate  the 
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RESIDUAL  COMPUTATION  SUMMARY 

ITERATION  NO.  MAX.  ELEMENT 

MAX.  TRACE 

0  0 . 9001E+00 

0 . 7547E-02 

1  0 . 3935E+00 

0 . 7807E+00 

2  0 . 1367E+00 

0 . 2728E+00 

3  0 . 1271E-01 

0 . 2540E-01 

4  0 . 1197E-03 

0 . 2389E-03 

5  0 . 7935E-06 

0 . 8964E-06 

EPSILON 

0 . 1000E-05 

MAXIMUM  NUMBER  OF  ITERATIONS 

NUMBER  OF  ITERATIONS  TO 

11 

CONVERGE  ON  EPSILON 

5 

CHOOSE  AN  INTEGER 

1.  CHANGE  EPSILON 

2.  CHANGE  MAX  ITERATIONS 

3.  RECOMPUTE  SPECTRAL  FACTOR 

4.  CONTINUE  PROGRAM 

Figure  3.2:  Residual  computation  summary  from  MIMO  case  benchmark  problem. 


E.  t ,  Es  (  l ,  1 ) 


Figure  3.3:  An  example  of  a  saturation  cost  vs.  tracking  cost  plot, 
design  frequency  response  data  with  this  k. 

The  final  objective  of  the  software  was  to  provide  a  means  for  displaying  and  obtaining 
output  of  the  frequency  response  data.  Specific  data  relevant  to  control  design  were  grouped 
into  a  plotting  menu  that  appears  in  figure  3.4.  These  data  can  be  plotted  on  the  screen 
or  as  a  hard  copy  in  the  various  configurations  shown  in  the  plot  type  menu  in  figure  3.5. 
In  addition,  there  is  a  feature  in  the  main  menu  ,that  permits  the  user  to  save  the  data  of 
figure  2.12  in  a  file  for  later  plotting. 

3.3.2  Testing  and  Benchmarks 

The  software  was  tested  using  four  benchmark  problems  -  a  pinned-  pinned  beam  with  a 
torque  at  one  end,  a  rigid  stick  on  a  cart  subject  to  a  disturbance,  and  a  flexible  stick  on 
a  cart  subject  to  a  disturbance  in  both  SISO  and  MIMO  cases.  These  benchmarks  will  be 
discussed  in  more  detail  in  a  following  section,  but  a  brief  summary  of  the  computational 
and  storage  requirements  for  these  problems  will  be  presented  here. 

The  software  was  tested  on  a  Micro  VAX  II,  and  the  CPU  time  for  each  benchmark 
problem  appears  in  Table  1.  This  time  represents  the  actual  CPU  time  required  to  perform 
all  of  the  necessary  computations  required  to  calculate  the  frequency  response  of  the  optimum 
controller.  As  can  be  seen  from  the  results,  there  was  a  substantial  increase  in  CPU  time  for 
the  MIMO  case  of  the  flexible  stick  on  a  cart,  which  is  a  reasonable  result  since  this  is  a  2 
by  2  matrix  case  whereas  the  other  benchmarks  were  scalar  cases. 

The  storage  requirements  for  each  benchmark  also  appear  in  Table  1.  The  frequency 
sampling  approach  is  very  data  intensive,  requiring  storage  of  many  arrays  containing  large 
numbers  of  elements.  Most  of  the  arrays  typically  contain  2048  elements  in  the  scalar  cases, 
and  8192  in  the  matrix  case. 

Additional  observations  of  the  results  of  the  benchmarks  show  that  accurate  results  of  the 
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PLOTTING  MENU. 


CHOOSE  A  NUMBER  BELOW. 


1.  XL(S) 

2.  YL(S) 

3.  OMEGA (S) 

4.  LAMBDA (S) 

5.  H(S) 

6.  HH(S) 

7.  Z(S) 


8.  K(S) 

9.  S(S),I-S(S) 

10.  L(S) 

11.  C(S) 

12.  ES.ET  COST  FUNCTIONS 

13.  PO(S) 

14.  RETURN  TO  MAIN  MENU 


Figure  3.4:  The  plotting  menu. 


PLOT  TYPE  MENU. 


CHOOSE  A  NUMBER  BELOW. 

1.  TYPE  OF  DATA  REPRESENTATION— -LINE 

2.  REAL  AND  IMAGINARY  VS.  FREQUENCY. 

3.  REAL  VS  IMAGINARY. 

4.  MAGNITUDE  AND  PHASE  VS  FREQUENCY. 

5.  RETURN  TO  PLOTTING  MENU. 


Figure  3.5:  The  plot  type  menu. 
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Benchmark 

CPU  Time 

Storage  Requirement 

Pinned-Pinned 

Beam 

6  min. 

2.7  Meg 

Rigid  Stick  on 
a  Cart 

5.5  min. 

2.7  Meg 

Flexible  Stick  on 
a  Cart  (SISO) 

9.5  min. 

2.7  Meg 

Flexible  Stick  on 
a  Cart  (MIMO) 

18  min. 

5.1  Meg 

Table  3.2:  Benchmark  computational  times  and  storage  requirements. 

frequency  response  curves  can  be  obtained  using  less  than  the  maximum  number  of  sample 
points  (1024)  which  reduces  the  computational  time  and  storage  requirements  substantailly. 
Figure  3.6  shows  the  controller  frequency  response  using  1024  points  and  figure  3.7  shows 
the  same  controller  using  only  256  points.  The  curves  are  very  similiar,  but  in  both  cases, 
the  low  frequency  end  of  the  respose  is  under-sampled  due  to  uniform  sampling.  However,  by 
reducing  the  bandwidth,  and  thereby  concentrating  more  of  the  sample  points  in  the  lower 
frequency  region  this  under-sampling  condition  can  be  reduced. 

3.3.3  Considerations  for  later  versions 

A  review  of  the  software  design  showed  that  a  substantial  decrease  in  memory  storage  can 
be  realized  by  a  reduction  in  the  number  and  size  of  temporary  storage  arrays.  Addition¬ 
ally,  the  size  of  some  of  the  permanent  arrays  can  also  be  reduced  by  rewriting  the  spectral 
factorization  code.  Furthermore,  the  program  can  broken  up  into  smaller  segments  that  per¬ 
form  intermediate  computations  and  then  store  pertinent  data  on  disk,  so  that  the  memory 
requirement  at  any  one  time  can  be  dramatically  reduced. This,  however,  would  increase  the 
software  run  time. 
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C(l,l) 


Figure  3.6:  Controller  frequency  response  with  1024  points. 


C(l.l) 


Figure  3.7:  Controller  frequency  response  with  256  points. 
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I 


control  torque 


y  -  point  of  observation 


Figure  3.8:  Pinned-Pinned  Beam  Control  Problem. 


3.4  Benchmark  Problems  Considered  in  Phase  1 
3.4.1  Control  of  a  pinned-pinned  beam. 

In  this  example  we  consider  a  control  system  design  problem  for  a  relativly  simple  structure 
with  well  known  transfer  function.  Consider  a  Bernoulli-Euler  beam  model  with  “pinned- 
pinned”  boundary  conditions  as  shown  in  the  figure.  The  beam  lateral  deformation  is  given 
by  y{t,z)  with  0  <  z  <  L  and  has  dynamics  described  by  the  PDE; 

M^-2(^i7^  +  K/0  =  O,  (3.44) 


with  boundary  conditions  at  z  =  0, 


(no  lateral  displacement) 


yM)  =  o, 


dz 2 


=  0, 


*=0 


(no  restraining  moment),  and  at  z  —  L, 


y{t,L)  =  0, 


d2y 

dz 2 


=  r, 


i=L 


where  the  control  moment  is  applied  at  the  right  hand  end  of  the  beam.  In  dimensionless 
form  the  PDE  can  be  written, 


V  _  2  >  d3y  d^y 

dt 2  ;  dtdz 2  dz 4 


(3.45) 


The  transfer  function  for  beam  control  is 

.  2  sin  sinh  X2 1  —  sin  X\  £  sinh  A2 


P(s,z)  =  L 


(AJ  +  A2)sin(Ai)sinh(A2) 


(3.46) 


SEI-89-03-15-WB 


A?  =  (-C  +  is/l^)sL\ 

=  «  +  V1  -  c2)^2> 


33 


where 


(3.47) 

(3.48) 


L  is  the  beam  length,  (  is  the  damping  factor,  and  0  <  z  <  L  is  the  observation  point  on 
the  beam.  The  resulting  transfer  function  is  meromorphic  with  poles  occurring  as 

P„  =  nVL2(-C  ±  iy/l-P)  (3-49) 


where  n  =  ±1,  ±2, .  .  ..  The  transfer  function  represents  a  stable  system  with  uniform  damp¬ 
ing  rate  given  by  (  and  P  6  A™- 

Numerically  stable  evaluation  of  the  beam  frequency  response  can  best  be  obtained  by 
evaluation  of  the  transfer  function  in  the  form 


_ .  x  L 2  /sinhA2f 

S'  =  (A?  +  A1)  V  sinh  A2 


sin 

sin  Ai 


(3.50) 


which  separates  the  exponential  terms  from  the  cyclic  terms  with  respect  to  multiplication. 
In  this  form  wide  band  frequency  response  data  can  be  obtained  with  sufficient  numerical 
precision. 

The  control  design  problem  we  examined  is  defined  by  a  choice  of  the  output  point  for 
regulation,  |  =  0.7,  the  length  of  the  beam  L  —  10.,  and  the  effective  damping  ratio, 
(  =  0.01.  We  note  that  since  P  is  pseudo-meromorphic  and  stable  the  coprime  factorization 
step  is  unnecessary.  We  take 


Nr  =  Nt  =  P,  Dr  =  Dl=  1. 

The  resulting  frequency  response  data  for  the  plant  model  is  shown  in  Figure  3.9.  Clearly, 
the  frequency  response  is  irrational  and  no  obvious  rational  approximation  is  evident. 

The  exogenous  inputs  describing  the  control  problem  are  given  as  follows:  the  output 
(load)  disturbance  is  assumed  constant  and  has  PSD  given  by  Gd{$)  =  the  feedback 
sensor  has  negligible  dynamics  (i.e.,  F0  =  F  =  1)  but  is  subject  to  noise  measurements 
characterized  by  PSD 

=  s4  +  2u>2(  1  —  2<,'2  ).s2  +  u>4 

The  control  problem  considered  is  output  regulation  of  the  beam  displacement  at  the  specified 
location  and  we  therefore  characterize  the  PSD  of  the  set  point  as  Ctu(^)  =  0.  The  PSD’s 
for  the  exogenous  inputs  are  shown  in  Figure  3.10. 

Control  design  procedes  by  evaluation  of  the  tradeoff  between  optimal  tracking  cost 
subject  to  a  constraint  on  the  saturation  of  some  critical  input  to  the  plant.  Here  we 
take  Q  =  1  so  that  saturation  cost  component  is  computed  with  respect  the  beam  control 
torque.  The  effective  optimal  cost,  J  =  Jt  +  kjsy  tradeoff  is  displayed  as  a  curve  of  tracking 
cost  Jt  vs.  saturation  cost  J$  in  Figure  3.13.  The  cost  components  is  estimated  from 
the  sampled  frequency  response  data  by  numerical  approximation  of  (2.57)-(2.58)  using 
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Figure  3.9:  Frequency  response  for  pinned-pinned  beam  control. 


Fraqufcf'icy 


Figure  3.10:  PSD  for  disturbance  and  sensor  noise  inputs  for  beam  control. 
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rectangular  quadrature.  The  sampled  data  was  computed  using  1024  uniform  sampling 
of  the  frequency  response  with  bandwidth  wbw  =  10.  The  performance  of  the  spectral 
factorization  steps  is  displayed  in  Figure  3.11.  The  algorithm  convergence  is  displayed  by 
examining  the  size  of  the  residual.  The  frequency  sampled  spectral  factors  for  the 

filtering  problem  and  A(s)  for  the  control  problem  are  displayed  in  Figure  3.12. 

To  check  the  sensitivity  of  the  computations  to  the  choice  of  bandwidth  we  recomputed 
the  optimal  control  tradeoff  study  using  1024  point  sampling  for  a  bandwidth  of  100.  The 
beam  frequency  response  is  shown  in  Figure  3.14.  The  tracking  vs.  saturation  cost  tradeoff 
for  this  bandwidth  is  displayed  in  Figure  3.15. 

The  tradeoff  in  control  system  performance  for  the  choice  of  tracking  vs.  saturation 
weighting  can  also  be  displayed  directly  in  the  frequency  domain  by  displaying  the  magnitude 
of  the  sensitivity  function  S(s)  and  its  complement  I  —  S(s).  Figure  3.16  displays  these 
functions  for  several  choices  of  the  parameter  k.  For  SISO  problems,  the  traditional  notions  of 
gain  and  phase  margin  can  be  readily  obtained  from  a  Nyquist  plot  of  the  loop  transmission, 
P(  ju>)C(  jut)  which  is  shown  for  the  1024  point  sampling  for  a  bandwidth  of  100  in  Figure  3.17. 

Once  the  tradeoff  analysis  of  saturation  vs.  tracking  is  resolved  the  optimal  controller 
is  specified  in  terms  of  its  frequency  response.  The  choice  of  k  =  0.1  acheives  a  saturation 
cost  for  this  problem  of  J9  —  7  representing  the  total  energy  in  the  actuating  torque  signal. 
The  resulting  controller  frequency  response  is  shown  in  Figure  3.18  for  1024  point  sampling 
of  over  a  bandwidth  of  100. 

All  frequency  response  plots  are  obtained  as  hardcopy  from  the  prototype  code  designed 
for  the  Phase  1  effort.  In  a  realistic  design  environment  the  frequency  response  data  can  be 
manipulated  and  displayed  interactively  to  support  tradeoff  studies  supporting  the  controller 
computational  phase  of  design.  In  Phase  2  we  propose  to  port  this  code  to  a  dedicated 
workstation  and  to  enhance  the  code  with  specific  provisions  for  tradeoff  analysis  of  various 
control  law  implementations. 

3.4.2  Vibration  Isolation  for  a  Simple  Elastic  Structure. 

A  prototype  MIMO  control  problem  for  active  vibration  isolation  is  considered  with  elastic 
dynamics  of  a  flexible  structure.  To  demonstrate  the  compatibility  with  available  models 
of  complex  flexible  structures  we  develop  a  simple  finite  element  representation  of  the  elas¬ 
tic  structure  response  and  use  the  resulting  finite  dimensional  model  to  approximate  the 
irrational  frequency  response. 

The  model  is  shown  in  Figure  3.19.  For  this  simple  problem  we  consider  all  motion 

constrained  in  the  x - z  plane.  The  model  consists  of  a  flexible  appendage  attached  via  a 

one  degree  of  freedom  rotational  joint  to  an  relatively  massive  carriage  which  is  supported 
mechanically  from  an  inertial  reference.  The  carriage  support  is  modeled  by  a  lumped 
stiffness  and  damping.  The  carriage  is  subject  to  a  disturbance  force  /.  Active  vibration 
isolation  is  achieved  in  this  simple  model  by  a  torque  r b  applied  at  the  rotational  joint 
between  the  flexible  appendage  and  the  carriage.  Additionally,  there  is  a  torque  ti  which  can 
be  applied  at  the  end  of  the  appendage.  Feedback  control  is  possible  using  measurements  of 
appendage  angle  0  and  angular  deformation  (shear  strain)  4>(L)  at  the  end  of  the  appendage. 
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RESIDUAL  COMPUTATION  SUMMARY 


ITERATION  NO. 
0 


MAX.  ELEMENT 
0.6426E-02 
0.1228E-04 
0.2389E-06 

OF  ITERATIONS 
TO 


MAX.  TRACE 
0 . 6426E-02 
0.1228E-04 
0.2389E-06 

0.1000E-05 

11 


2 

EPSILON 

MAXIMUM  NUMBER 
NUMBER  OF  ITERATIONS 
CONVERGE  ON  EPSILON 


Figure  3.11:  Convergence  of  filter  spectral  factorization  for  pinned- pinned  beam. 
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Figure  3.12:  Spectral  Factors  for  pinned-pinned  beam  control. 
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E  L  ,  Es (  1  .  1 ) 


Figure  3.13:  Tracking  vs.  Saturation  Tradeoff  for  Beam  Control. 


P(i.i) 


Figure  3.14:  Frequency  response  for  pinned-pinned  beam  control. 


SEI-89-03-15-WB 


Figure  3.15:  Tracking  vs.  Saturation  Tradeoff  for  Beam  Control. 
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Figure  3.16:  Design  Tradeoff  using  Sensitivity  Function  for  Beam  Control. 
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Figure  3.17:  Nyquist  Plot  of  Beam  Control  Loop  for  k  =  .1. 


Figure  3.18:  Optimal  Controller  Frequency  Response  for  k  =  .1. 
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model  parameter 

value 

description 

mc 

0.9 

mass  of  the  carriage 

k 

1.0 

support  stiffness 

b 

1.0 

support  damping 

L 

15. 

appendage  length 

P 

0.0707 

mass  density  of  appendage 

A 

0.0942 

area  cross  section  for  appendage 

la 

5.302 

area  moment  of  inertia  for  appendage 

E 

16.0 

modulus  elasticity  of  appendage 

kG 

6.4 

effective  shear  modulus 

Ci 

C2 

0.01206 

1.697 

damping  constants  for  material  dissipation 

Table  3.3:  Model  Parameters  for  Vibration  Isolation  Problem. 


Thus  the  plant  transfer  function  P(s)  is  2  x  2. 

The  dynamic  model  is  obtained  by  application  of  Hamilton’s  principle  in  terms  of  the  gen¬ 
eralized  coordinates  {z,  0,  r/(z),  <f>(z)}  where  x  is  the  horizontal  displacement  of  the  carriage, 
9  is  the  angular  displacement  of  the  appendage  at  the  joint,  r/(z)  is  the  lateral  deformation 
of  the  appendage  relative  to  the  centerline  given  by  9,  and  4>{z)  is  the  angular  deformation 
of  the  appendage  cross  section  over  the  length  of  the  appendage;  0  <  z  <  L.  The  model 
parameters  are  given  in  the  Table  3.3.  Parameters  are  dimensionless  and  chosen  to  provide 
reasonably  well  scaled  spectral  response  of  typical  flexible  structure  response. 

Under  the  above  assumptions  the  system  kinetic  energy  takes  the  form, 

T(x,  0,  r), <t>)  =  ^mcx2  +  ^J  {M(i  +  i{z)  -  z9)2  +  pla{9  +  <j>{z)f}dz ,  (3.51) 


and  the  system  potential  energy  is, 


'dj 

dz 


V(«,  9,VA)  =  +  \I0L\EI‘(%)  +*0-4 

We  also  assume  a  dissipation  function  of  the  form, 


—  (f> )  >  dz. 


(3.52) 


(3.53) 


which  models  material  losses  which  are  expected  to  be  dominant  in  terms  of  dissipation  for 
space  applications. 

A  finite  dimensional  model  can  be  obtained  by  the  Finite  Element  Method  (FEM)  by 
introducing,  for  example,  assumed  modes  and  expansions  of  the  form, 

r/(t,  z )  «  $r(z)7/(t),  <f>(t ,  z)  %  $T(z)4>(t), 

where  t},  <f>  are  N- vectors  of  FEM  degrees  of  freedom.  We  prefer  the  introduction  of  “assumed 
modes”  from  FEM  using  collocation  by  splines.  We  note  that  under  the  above  assumptions 
1st  order  (linear)  splines  are  sufficient  [28]. 
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Applying  the  above  FEM  approximation  to  the  energy  functions  and  solving  the  Euler- 
Lagrange  equations  for  the  reduced  system  Lagranian  L  =  T  —  V  obtains  the  dynamic 
equations  in  the  form, 

Mw(t)  +  Bw(t)  +  Iiw(t)  =  Eu(t),  (3.54) 

y(t)  =  Ctw(t)  + 

where  w  —  [x,  9,  ij,  <j>]T ,  u  =  [■q,,r^,/]T,  and  y  =  [9,(f>(L)]T.  (Model  data  in  this  form  are 
typically  obtained  from  a  variety  of  standard  FEM  codes  such  as  NASTRAN.)  The  matrix 
parameters  have  the  form 


M  = 


— 

Nex 

nL 

T)X 

0 

-n9x 

Je 

N9<p 

Nnx 

— 

N9n 

Nr, 

0 

0 

Ne<p 

0 

N't, 

'  b 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Br, 

Br rt> 

1 

0 

0 

B^ 

Bj, 

'  k 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Kr, 

Kr,<j> 

0 

0 

K^ 

iu 

(3.55) 


(3.56) 


(3.57) 


Here  Mtot  is  the  system  total  mass,  and  J$  is  the  effective  rigid  inertial  moment  for  the  0  mo¬ 
tion.  The  other  model  parameters  arise  from  the  FEM  approximation  and  are  characterized 
by,  for  example,  scalar  expressions  of  the  form 


Nex 


z  dz , 


row  matrices  of  the  form, 

NTn,  =  f\AiT(z)Jz, 
and  tridiagonal  matrices  of  the  form, 


Nr, 

Kr, 


fL  pAi(z)^T(z)dz, 
Jo 

fL  d$d$T  , 

/  kG~z — dz 
Jo  OZ  oz 


(cf  [28]  for  details  on  the  development  of  FEM  models  by  collocation  by  splines.)  The  matrix 
E  is  (2  +  2 N)  x  3  with  £7(1,1)  =  1,  £7(2,3)  =  1,  and  £7(2  +  2N,2)  =  1.  The  composite 
matrix  [Ci,  Cj]  is  the  transpose  of  the  first  two  columns  of  £7. 
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The  parameters  of  the  model  were  chosen  to  provide  a  spectrum  for  the  numerical  com¬ 
putations  of  the  frequency  response  characteristic  of  flexible  space  structures  with  damping 
arising  primarily  from  material  losses  and  to  illustrate  opportunities  for  tradeoff  in  the  com¬ 
putation  of  stable  coprime  factorization.  We  have  also  structured  the  control  problem  by 
choosing  collocated  sensing  and  actuation.  Figure  3.20  illustrates  the  position  of  the  poles 
and  transmission  zeros  for  the  2x2  transfer  function  P(s)  which  takes  inputs  (r^rx,)7^  to 
outputs  y.  The  effect  of  the  disturbance  loading  is  given  by  the  transfer  function  P0(s)  from 
/  to  y. 

We  note  that  the  system  model  is  unstable  with  one  real  pole  in  the  right  half  plane. 
The  unstable  mode  is  associated  with  the  unstable  reaction  of  the  flexible  appendage  which 
is  attached  to  the  carriage  by  a  one  degree  of  freedom  rotational  joint  with  no  dynamic 
restraint.  It  will  be  required  to  obtain  a  stable  coprime  factorization  for  the  system  trans¬ 
fer  function  before  proceding  with  the  numerical  frequency  sampling  computations  of  the 
controller  frequency  response. 

Stable  coprime  factorization  can  be  directly  obtained  from  the  model  (3.54),  or  equiva¬ 
lently,  by  transforming  the  model  to  state  space  form, 

x  =  Ax  +  Bu ,  (3.58) 

y  =  Cx, 


with 


0  I 

0 

—M~lK  —M~1B 

II 

2Q 

M~lE 

c  =  [CuC2\. 


(3.59) 


As  discussed  in  the  previous  section  stable  coprime  factorizations  can  be  obtained  from 
the  state  space  realization  (cf.  (3.2)-(3.4)).  First,  find  real  matrices  H  and  F  so  that  the 
eigenvalues  of  Ah  =  A  +  HC  (resp.  A?  =  A  +  BF )  are  in  the  left  half  plane.  Then  the 
coprime  factors  can  be  given  as, 

[Nt,  Dt)(s)  =  [0,  Ip]  +  C[sl  -  Ah]~1[B,  H],  (3.60) 

{sI-Af}-'B.  (3.61) 

The  pole/zero  plots  for  the  respective  factors  are  shown  in  Figure  3.21  where  we  have  chosen 
the  poles  of  the  factors  by  translating  the  real  part  of  the  unstable  and  nearly  unstable  poles 
(including  the  real  pole  at  s  =  1.73  and  the  double  pole  at  the  origin)  of  P(s)  to  the  left. 
From  these  figures  it  is  clear  that  the  transmission  zeros  of  N(  (resp.  D()  are  the  transmission 
zeros  (resp.  poles)  of  P.  The  poles  of  [Afy,  Di ]  represent  the  stabilizing  choices  of  the  open 
loop  poles. 

The  disturbance  force  is  modeled  via  a  narrowband  PSD  centered  about  cu  =  1  given  by 

_<,2 

Gf^  ~  7*  +  2(1  -  2 c2)s2  +  r 


Nr 

Dr 


+ 


C 
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Figure  3.19:  Vibration  Isolation  Model  with  Flexible  Structure. 
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Figure  3.20:  Pole/Zero  plot  for  Vibration  Isolation  Model. 
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Figure  3.21:  Pole/zero  Locations  for  the  Coprime  Factors. 

This  might  represent  a  periodic  motion  of  some  onboard  subsystem  components  for  a  space 
craft  example.  Assume  the  displacement  sensors  are  subject  to  independent,  zero  mean  noise 
processes  described  by  PSD’s  of  the  form 

54  +  20,2(1  -2C)s2 +  U,*' 

Tradeoff  analysis  of  tracking  vs.  saturation  was  performed  and  a  value  of  k  =  .5  was  cho¬ 
sen  for  the  candidate  design.  The  resulting  closed  loop  control  bandwidth  can  be  determined 
from  the  frequency  response  of  the  2x2  sensitivity  function  5  and  its  complement  I  —  5. 
In  Figure  3.22  the  singular  values  of  these  matrix  frequency  responses  are  shown  plotted  as 
dB  gain  vs.  log  frequency.  From  these  figures  one  can  obtain  an  estimate  of  the  desired 
closed  loop  control  bandwidth  obtained  by  Wiener-Hopf  optimization.  Note  the  numerical 
noise  evident  in  the  singular  value  plots  of  I  —  S  for  high  frequencies.  This  is  a  result  of 
the  single  precision  computations  currently  implemented  in  the  software.  We  display  the 
specification  of  the  optimal  control  law  in  terms  of  the  computed  frequency  samples  of  its 
scalar  components  via  the  Bode  plots  in  Figure  3.23. 

We  remark  that  the  problem  considered  here  is  characteristic  of  structural  vibration 
control  problems  in  several  ways.  One  important  feature  is  evident  by  examination  of  the 
frequency  response  of  the  open  loop  plant  transfer  function.  In  Figure  3.24  we  plot  the 
minimum  and  maximum  singular  values  of  the  frequency  response  of  P,  The  rigid  body  mode 
at  the  origin  is  evident  in  the  effective  gain  (i.e.  maximum  singular  value)  increasing  without 
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bound  for  decreasing  frequency.  However,  the  ratio  of  maximum  to  minimum  singular  values 
is  also  increasing  without  bound  for  decreasing  frequencies.  This  indicates  a  numerical 
singularity  of  the  effective  frequency  response  for  5  =  0.  But  this  is  characteristic  of  the 
redundant  actuators  and  sensors  for  the  effective  rigid  (i.e.  low  frequency)  response  of  the 
appendage.  It  is  interesting  to  note  that  even  with  such  ill-conditioned  transfer  function 
model  that  the  numerical  algorithms  for  spectral  factorization,  causal  projection,  and  other 
components  of  the  Wiener-Hopf  control  design  computations  produce  numerically  stable 
results.  We  intend  to  direct  software  and  algorithm  development  in  the  Phase  2  proposal 
toward  the  development  of  various  options  for  extended  precision  computations  and  options 
for  MIMO  control  loop  shaping  for  transfer  function  models  with  low  frequency  redundancies 
(i.e.  near  singular  numerical  behavior). 

3.5  Realtime  Control  Law  Implementation 

The  frequency  sampled  computations  obtain  a  specification  for  the  frequency  response  of 
the  ideal  (optimal)  controller  via  its  sampled  representation.  The  design  engineer  now  has 
several  options  for  implementing  the  control  law.  Typically  for  aerospace  applications,  high 
speed  digital  computers  will  be  needed  for  real  time  control. 

A  general  representation  of  the  dynamic  control  law  is  in  terms  of  a  convolution, 

u{t)  =  c{t)  *  y(t), 

of  the  measurements  y(t)  with  a  controller  impulse  response,  c(t).  We  contend  that  the 
technology  for  implementing  high  speed  convolution  or  filtering  for  real  time  control  of  flexible 
structures  is  now  available  and  can  be  realized  using  digital  signal  processing  methods. 

The  convolution  above  can  be  approximated  to  any  desired  degree  using  a  discrete  time 
representation, 

u(nT)  =  c((n  -  k)T)y(nT). 

k= 0 

In  general  such  a  system  can  be  represented  as  a  causal  recursive  linear  system ,  i.e. ,  one 
whose  present  output  is  computed  from  past  and  present  inputs  and  outputs.  However,  in 
many  cases  the  response  can  be  approximated  by  a  nonrecursive  system ;  i.e.,  one  whose 
present  output  depends  only  on  past  and  present  inputs.  Thus  a  nonrecursive  realization  (if 
one  exists)  would  take  the  form  of  a  moving  average, 

m 

u{nT)  =  Y2  aky((n  -  k)T). 
k= o 

Such  a  realization  is  often  called  a  Finite  Impulse  Response  (FIR)  filter,  since  the  effective 
impulse  response  will  become  zero  after  m  time  steps.  If  the  ideal  continuous  controller  has 
transfer  function  which  is  analytic  in  the  closed  right  half  plane  and  is  strictly  proper  then  its 
impulse  response  approachs  zero  as  t  — >  oo.  Such  an  impulse  response  can  be  approximated 
arbitrarily  well  by  a  FIR  filter. 
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Figure  3.22:  Singular  Value  plots  of  the  Sensitivity  Function  for  k  =  .5. 
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Figure  3.23:  Bode  plots  for  the  scalar  elements  of  C(s) 
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Figure  3.24:  Singular  value  plots  of  P(s). 

Frequency  Sampling  Filter  realization.  An  example  of  an  FIR  filter  is  the  frequency 
sampling  filter  whose  frequency  response  is  chosen  to  interpolate  a  desired  frequency  response 
by  exactly  matching  a  set  of  N  frequency  samples.  Consider  the  case  where  the  frequency 
samples  are  uniformly  spaced.  Since  the  filter  is  implemented  in  discrete  time  its  transfer 
function  can  be  described  using  z-transforms.  The  frequency  response  samples  are  therefore 
given  as; 

zk  =  ejk2,t/N, 

for  k  =  0, . . . ,  N  —  1,  (the  N  roots  of  unity)  which  correspond  to  the  frequencies 

uik  =  ku),/N, 

where  u;,  is  the  time  sampling  rate  for  the  discrete  time  realization. 

The  interpolating  function  in  the  z-plane  can  be  given  as 

1  -  z~N 

=  JV(1  -ejk2w/Nz-l)' 

with  k  =  0, . . . ,  N  —  1.  This  function  has  the  interpolating  property; 

F,  ( Jt2n/N\  —  /  0  f°r  l  7^  & 

fcl  ’  ~  \  1  for  l  =  k 

If  we  take  the  desired  frequency  samples  as  Ck  then  the  transfer  function 

c(>)=Ew,(.) 

fc=o 
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will  realize  the  frequency  sampling  filter.  It  is  easy  to  show  that  the  frequency  sampling 
filter  has  transfer  function 

1  _  r~N  7V_1  n  N~l 

G(z)  =  ~~N~  ^ 

iy  k=o  1  c  z  *=o 

where  t  he  coefficients  are  given  as 

c,=  jE 

k=0 

for  (  =  0,...,  N  -  1. 

The  coefficients  c0, . . . ,  c^-\  are  known  as  the  Inverse  Discrete  Fourier  Transform  (IDFT) 
of  the  sequence  Co, . . . ,  CV-i-  Thus  a  simple  implemenation  of  the  frequency  sampling  filter 
which  has  been  used  in  a  wide  range  of  digital  signal  processing  applications  is  described  as 
follows.  Sample  the  input  time  waveform  y(t)  with  uniform  sampling  rate  and  sequencially 
load  the  samples  in  a  buffer.  After  every  k  <  N  samples  are  stored  in  the  input  buffer 
the  convolution  is  implemented  by  performing  a  DFT  of  the  input  sequence.  Then  perform 
IDFT  of  the  product  of  this  result  with  the  frequency  samples  (i.e.,  the  coefficients  Ck) 
to  obtain  the  sampled  time  representation  of  the  output  y(t ).  The  output  buffer  is  then 
sequentially  clocked  to  the  output  D/A  channel.  This  illustrates  that  the  delay  in  processing 
is  dependent  on  the  computational  delay  for  DFT  processing  of  N  length  sequences.  It 
is  well  known  that  the  convolution  operation  just  described  is  most  efficiently  implemented 
using  the  Fast  Fourier  Transform  (FFT) — a  special  form  of  the  DFT  optimized  for  minimum 
number  of  multiplications  [26]. 

In  general,  the  problem  of  realizing  closed  loop  controller  by  FIR  filters  is  complicated 
by  the  available  processing  speeds  and  the  duration  of  the  ideal  continuous  time  controller 
impulse  response.  The  construction  of  stablizing  controllers  via  stable  coprime  factorization 
can  potentially  play  a  significant  role  in  practical  aspects  of  FIR-based  controller  realization. 

Consider  a  general  stabilizing  feedback  controller  given  in  terms  of  the  relation  (2.41). 
The  control  law  can  be  expressed  in  the  frequency  domain  as 

ti(s)  =  C(s)y(s) 

or  equivalently, 

u(s)  =  [y(s)  +  Dr(s)K(s)}C(s)  (3.62) 

A'(s)£(5)  =  y(s)  +  Nr{s)K(s)t(s),  (3.63) 

where  we  take  A",  N ,  D,  K  can  be  chosen  to  be  “arbitrarily  stable^.  By  this  we  mean  that 
their  respective  poles  can  be  chosen  to  be  to  the  left  of  a  vertical  axis  in  the  complex  plane 
whose  real  part  can  be  chosen  arbitrarily  far  into  the  left  half  plane.  The  controller  (7(<s) 
may  not  be  arbitrarily  stable  however  since  X~l  may  have  have  poles  to  the  right  of  the 
arbitrarily  chosen  verticle  line.  Indeed,  in  some  cases  C  may  not  be  stable  at  all  in  which 
case  implementation  of  the  control  law  by  direct  FIR  processing  is  not  feasible.  However,  it 
is  easy  to  show  that  realization  can  be  obtained  by  feedback. 
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Figure  3.25:  Option  for  Stable  FIR  controller  implementation 

To  see  how  this  can  be  we  obtain  X  in  the  form, 

X  =  U~  X 

where  U  has  an  “stable”  inverse  (i.e.  has  all  poles  to  the  left  of  the  arbitrary  verticle  line) 
and  A'  has  only  “unstable”  poles.  Note  that  for  most  flexible  structure  models  A’  will  be 
rational.  Then  we  can  write  the  control  law  in  the  form, 

0  \Y  +  D,K\  ](v\_r(>\ 

U-'  {U-'N.K  -  /]  J  \  f  )  \(  I 

which  can  be  implemented  as  a  stable  filter  with  external  feedbacks  as  shown  in  Figure  3.25. 
Note  that  each  component  of  this  feedback  realization  of  the  control  law  is  arbitrarily  stable. 

We  have  shown  that  the  controller  can  be  realized  by  a  feedback  connection  of  compo¬ 
nents  each  of  which  has  impulse  response  which  can  be  tailored  for  specific  implementation 
requirements  of  the  components  available  to  the  designer.  In  particular,  each  component 
is  stable  and  can  be  realized  by  an  FIR  filter.  In  Phase  2  we  will  focus  on  implemenation 
considerations  for  FIR  type  processing  and  the  options  for  feedback  realizations  of  this  type 
seems  to  offer  increased  flexibility  for  such  designs. 

4  Conclusions  and  Directions 

Results  from  the  Phase  1  study  have  demonstrated  that  the  computational  problem  of  fre¬ 
quency  domain  design  can  be  solved  using  frequency  sampling.  As  such  the  computations 
can  be  readily  extended  to  certain  classes  of  irrational  transfer  functions;  such  as  typically 
arise  in  control  of  flexible  structures.  Advantages  of  the  approach  include: 

1.  The  computation  of  control  by  frequency  response  sampling  can  embrace  design  prob¬ 
lems  consisting  of  both  irrational  and  rational  transfer  function  models  using  the  same 
algorithms. 

2.  The  required  frequency  response  data  for  model  building  can  be  obtained  from  ana¬ 
lytical  models  of  the  elastic  dynamics,  finite  element  models  of  structural  responses, 
and/or  empirical  data  obtained  from  measurements  on  full  or  partial  scale  models. 

3.  The  flexiblitity  of  control  law  implementation  is  left  to  the  design  engineer. 
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We  believe  that  a  natural  approach  for  flexible  structure  control  law  implementation  is  to 
realize  the  required  processing  (i.e.,  real  time  convolution)  by  discrete  time  methods  and  to 
implement  using  advanced  high  speed  VLSI  technology.  For  real  time  control  applications 
the  processing  delay  incurred  in  FIR  implementation  is  a  critical  factor  which  may  effect 
the  system  stability  margins.  The  major  goal  of  the  phase  2  effort  will  be  to  demonstrate 
that  high  speed  implementations  can  reduce  the  processing  delay  to  acceptable  levels  for 
application  to  control  of  flexible  structures. 

Another  important  issue  in  realtime  implementation  arises  when  the  ideal  controller 
(7(s)  is  open  loop  unstable  (i.e.,  has  poles  in  the  closed  right  half  plane).  In  such  cases  the 
impulse  response  can  not  be  directly  approximated  by  an  FIR  filter.  However,  we  have  shown 
that  the  impulse  response  can  be  realized  by  a  feedback  interconnection  of  FIR  filters  [18]. 
Moreover,  the  flexibility  inherent  in  the  algebraic  description  of  all  stabilizing  controllers  for 
a  given  plant  transfer  function  suggests  several  new  options  for  controller  realization  by  FIR 
processing  which  we  will  propose  to  investigate  further  in  Phase  2. 
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